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Summary 


This research program deals with the application of high-performance computing methods for the 
analysis of complete jet engines. We have initited this program by applying the two-dimensional 
parallel aeroelastic codes to the interior gas flow problem of a by-pass jet engine. The fluid mesh 
generation, domain decomposition and solution capabilities were successfully tested. We have then 
focused attention on methodology for the partitioned analysis of the interaction of the gas flow with 
a flexible structure and with the fluid mesh motion that results from these structural displacements. 
This is treated by a new ALE technique that models the fluid mesh motion as that of a fictitious mass- 
spring network. New partitioned analysis procedures to treat'this coupled 3-component problem 
are developed. These procedures involved delayed corrections and subcycling. Preliminary results 
on the stability, accuracy and MPP computational efficiency are reported. 






1. OVERVIEW 


The present program deals with the application of high-performance parallel computation for the 
analysis of complete jet engines, considering the interaction of fluid, thermal and mechanical 
components. The research is driven by the simulation of advanced aircraft propulsion systems, 
which is a problem of primary interest to NASA Lewis. 

The coupled problem involves interaction of structures with gas dynamics, heat conduction and heat 
transfer in aircraft engines. The methodology issues to be addressed include: consistent discrete 
formulation of coupled problems with emphasis on coupling phenomena; effect of partitioning 
strategies, augmentation and temporal solution procedures; sensitivity of response to problem 
parameters; and methods for interfacing multiscale discretizations. The computer implementation 
issues to be addressed include: parallel treatment of coupled systems; domain decomposition and 
mesh partitioning strategies; data representation in object-oriented form and mapping to hardware 
driven representation, and tradeoff studies between partitioning schemes with differing degree of 
coupling. 


2. STAFF 

Two graduate students have been partly supported by this grant. M. Ronaghi (U.S. citizen) began his 
graduate studies at Colorado on January 1993. He has completed a M.Sc. in Aerospace Engineering 
on May 1994 but plans to pursue his Ph. D. in Mechanical Engineering under a different research 
program that fits his background and experience better. 

U. Gumaste (permanent U.S. resident) began his graduate studies at Colorado in the Fall semester, 
but worked in this project during June and July 1993 as an hourly research assistant. Mr. Gumaste 
has a B.Tech in Civil Engineering from the Indian Institute of Technology, Bombay, India. He 
completed his Ph. D. course requirement this semester with the transfer of graduate credit units 
from the University of Maryland. During this period he was partly supported as a Teaching Assistant. 
He plans to familiarize himself with our aeroelastic codes during May-June and will visit NASA 
Lewis for 6 weeks during July and August. 

The methods development for this project has been greatly benefitted from the presence of three 
visiting Post-Docs. Dr. S. Lanteri conducted extensive experimentation on several computational 
algorithms for compressive viscous flow simulation on the iPSC-860, CM-5 and KSR- 1 as reported 
in Appendix II. Dr. N. Maman has implemented “mesh matching” techniques that connect separately 
generated fluid and structural meshes. Dr. S. Pipemo has developed and evaluated implicit and 
subcycled partitioned analysis procedures for the interaction of structure, fluid and fluid-mesh 
motion. A new approach to augmentation of the governing semi-discrete equations that improves 
stability while keeping communications overhead modest was investigated. Results from this study 
are presented in Appendix I. 
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3. DEVELOPMENT OF PARTITIONED ANALYSIS METHODS 


The first parallel computations of a jet engine, presented in the first Progress Report, dealt with 
the fluid flow within a jet engine structure that was considered rigid and hence provides only 
guiding boundary conditions for the gas flow. When the structural flexibility is accounted for two 
complications occur: 

1. The engine simulation algorithm must account for the structural flexibility though periodic 
transfer of interaction information, and 

2. The fluid mesh must smoothly follow the relative structural motions through an ALE (Adaptive 
Lagrangian Eulerian) scheme. The particular ALE scheme selected for the present work makes 
use of Batina’s proposed pseudo-mechanical model of springs and masses overlaid over the 
fluid mesh. 

Research work during the period July 1993 through January 1994 was dominated by the treatment 
of two subjects: partitioned analysis of fluid-structure interaction (FSI) and accounting for fluid 
mesh motions. The partitioned analysis algorithm developed for the FSI problem is always implicit 
in the structure (because of its larger time scale of significant vibratory motions) and either explicit 
or implicit for the gas flow modeled by the Navier-Stokes equations. Subcycling, in which the 
integration stepsize for the fluid may be smaller than that used in the structure, was also studied. 

General Requirements 

The fundamental practical considerations in the development of these methods are: (1) numeri- 
cal stability, (2) fidelity to physics, (3) accuracy, and (4) MPP efficiency. Numerical stability is 
fundamental in that an unstable method, no matter how efficient, is useless. There are additional 
considerations: 

1. Stability degradation with respect to that achievable for the uncoupled fields should be min- 
imized. For example, if the treatment is implicit-implicit (I-I) we would like to maintain 
unconditional stability. If the fluid is treated explicitly we would like to maintain the same 
CFL stability limit. 

2. Masking of physical instability should be avoided. This is important in that flutter or diver- 
gence phenomena should not be concealed by numerical dissipation. For this reasons all time 
integration algorithms considered in this work must exclude the use of artificial damping. 

Stability vs. Communication-Overhead Tradeoff 

The degradation of numerical stability degradation is primarily influenced by the nature of infor- 
mation exchanged every time step among the coupled subsystems during the course of partitioned 
integration. A methodology called augmentation that systematically exploits this idea was devel- 
oped by Park and Felippa in the late 1970s. The idea is to modify the governing equations of one 
subsystem with system information from connected subsystems. The idea proved highly successful 
for the sequential computers of the time. A fresh look must be taken to augmentation, however, in 
light of the communications overhead incurred in massively parallel processing. For the present 
application three possibilities were considered: 
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No augmentation. The 3 subsystems (fluid, structure and ALE mesh) exchange only minimal 
interaction state information such as pressures and surface-motion velocities, but no information 
on system characteristics such as mass or stiffness. The resulting algorithm has minimal MPP 
communication overhead but poor stability characteristics. In fact the stability of an I-I scheme 
becomes conditional and not too different from that of a less expensive I-E scheme. This degradation 
in turns significantly limits the stepsize for both fluid and structure. 

Full augmentation. This involves transmission of inverse-matrix-type data from one system to an- 
other. Such data are typified by terms such as a a structure -to-fluid coupling-matrix times the inverse 
of the structural mass. Stability degradation can be reduced or entirely eliminated; for example I-I 
unconditional stability may be maintained. But because the transmitted matrix combinations tend 
to be much less sparse than the original system matrices, the MPP communications overhead can 
become overwhelming, thus negating the benefits of improved stability characteristics. 

Partial augmentation. This new approach involves the transmission of coupling matrix information 
which does not involve inverses. It is efficiently implemented as a delayed correction to the 
integration algorithm by terms proportional to the squared stepsize. The MPP communication 
requirements are modest in comparison to the fully-augmented case, whereas stability degradation 
can be again eliminated with some additional care. 

The partial augmentation scheme was jointly developed by S. Pipemo and C. Farhat. Its derivation 
is fully reported in Appendix I, which has been submitted for publication. The general methodology 
is first applied to a staggered FSI algorithm (staggering is a special form of partition) with common 
timestep, and then extended to cover subcycling. The reduced communications overhead has been 
recently verified on simple problems with preliminary tests on a iPSC860 Hypercube. Tests of this 
new scheme for more complex geometries will be carried out this summer. 

Algorithmic Effects of Dynamic Fluid Mesh 

The first one-dimensional results on the effect of a dynamic fluid mesh on the stability and accuracy 
of the staggered integration were obtained by C. Farhat and S. Pipemo and are also discussed in 
Appendix I. A doctoral student, M. Lesoinne, supported by NSF is extending these calculations to 
the multidimensional case. 

Appendix II contains a report by Charbel Farhat and S. Lanteri that summarizes recent experiences 
with the application of the two-dimensional Navier Stokes equations combined with the pseudo- 
mechanical network approach to dynamic mesh motion. Performance results obtained on the iPSC- 
860, KSR- 1 and CM-5 parallel computers are reported for both fixed and moving fluid meshes. 

4. FUTURE WORK 

The present work undertaken since the renewal of the grant on March 1 994 is focused on axisym- 
metrizing the fluid-structure interaction codes, in ‘which the new partially-augmented staggered 
schemes and mesh motion techniques are implemented. We still need to introduce more physical 
effects in the gas flow, namely compression, diffusion and combustion. The modeling experience 
that Mr. Gumaste will acquire during his visit to NASA Lewis should help with the prosecution of 
these items. 
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January 1994 


Abstract 

In order to predict the dynamic response of a flexible structure in a fluid flow, the 
equations of motion of the structure and the fluid must be solved simultaneously. In this 
paper, we present several partitioned procedures for time-integrating this focus coupled 
problem and discuss their merits in terms of accuracy, stability, heterogeneous com- 
puting, I/O transfers, subcycling, and parallel processing. All theoretical results are 
derived for a one-dimensional piston model problem with a compressible flow, because 
the complete three-dimensional aeroelastic problem is difficult to analyze mathemati- 
cally. However, the insight gained from the analysis of the coupled piston problem and 
the conclusions drawn from its numerical investigation are confirmed with the numeri- 
cal simulation of the two-dimensional transient aeroelastic response of a flexible panel 
in a transonic nonlinear Euler flow regime. 


1. Introduction 

In order to predict the dynamic response of a flexible structure in a fluid 
flow, the equations of motion of the structure and the fluid must be solved si- 
multaneously. One difficulty in handling numerically the fluid/structure coupling 
stems from the fact that the structural equations are usually formulated with 
material (Lagrangian) co-ordinates, while the fluid equations are typically writ- 
ten using spatial (Eulerian) co-ordinates. Therefore, a straightforward approach 
to the solution of the coupled fluid/structu,re dynamic equations requires mov- 
ing at each time-step at least the portions of the fluid grid that are close to the 
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moving structure. This can be appropriate for small displacements of the struc- 
ture but may lead to severe grid distorsions when the structure undergoes large 
motion. Several different approaches have emerged as an alternative to partial 
re-gridding in transient aeroelastic computations, among which we note the ar- 
bitrary Lagrangian-Eulerian (ALE) formulation [1-3], the co- rotational approach 
[4,5], and dynamic meshes [6] (see also [28] for a review). All of these approaches 
treat a computational aeroelastic problem as a coupled two-field problem. 

However, the moving mesh itself can be formulated as a pseudo-structural 
system with its own dynamics [7], and therefore, the coupled transient aeroelastic 
problem can be formulated as a three- rather than two-field problem: the fluid, 
the structure, and the dynamic mesh. The semi-discrete equations governing this 
three-way coupled problem can be written as follows: 

^(A(x,t)W(x,t)) + F'(W(x,t),x,x) = F d (W(x,t)) 

= /“WM)) (1) 

~ d 2 x ~dx - 

M W i ' D dt + Kx = h< q 

where x is the position of a moving fluid grid point, W is the fluid state vector, 
A results from the finite element/volume discretization of the fluid equations, 
F£ is the convected vector of convective fluxes [7], F d is the vector of diffusive 
fluxes, q is the structural displacement vector, f tnt denotes the vector of internal 
forces in the structure, f ext the vector of external forces, M is the finite element 
mass matrix of the structure, M, D, and K are fictitious mass, damping, and 
stiffness matrices associated with the fluid moving grid and constructed to avoid 
any parasitic interaction between the fluid and its grid, or the structure and the 
moving fluid grid [7], and K c is a transfer matrix that describes the action of the 
motion of the structural side of the fluid/structure interface on the fluid dynamic 
mesh. For example, M = D = 0, and K = I corresponds to a rigid mesh motion 
of the fluid grid around an oscillating airfoil, and M = D = 0 corresponds to the 
spring-based mesh motion scheme introduced in [6]. 

Each of the three components of the coupled problem described by Eqs. (1) 
has different mathematical and numerical properties, and distinct software imple- 
mentation requirements. For Euler and Navier-Stokes flows, the fluid equations 
are nonlinear. The structural equations may be linear or nonlinear. The semi- 
discrete equations governing the pseudo-structural fluid grid system axe linear. 
The matrices resulting from a linearization procedure are in general symmetric 
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for the structural problem, but they are typically unsymmetric for the fluid prob- 
lem. Morevoer, the nature of the coupling in Eqs. (1) is implicit rather than 
explicit, even when the fluid mesh motion is ignored. The fluid and the structure 
interact only at their interface, via the pressure and the motion of the physical 
interface. However, the pressure variable cannot be easily isolated neither from 
the fluid equations nor from the fluid state vector W. Consequently, the numeri- 
cal solution of Eqs. (1) via a fully coupled monolithic scheme is computationally 
challenging and software-wise unmanageable. 

Alternatively, Eqs. (1) can be solved via partitioned procedures [8]. This 
approach offers several appealing features including the ability to use well estab- 
lished discretization and solution methods within each discipline, simplification 
of software development efforts, and preservation of software modularity. Tra- 
ditionally, transient aeroelastic problems have been solved via the simplest pos- 
sible partitioned procedure whose cycle can be described as follows: a) advance 
the structural system under a given pressure load, b) update the fluid mesh ac- 
cordingly, and c) advance the fluid system and compute a new pressure load 
[9-12]. Occasionally, some investigators have advocated the introduction of a few 
predictor-corrector iterations within each cycle of this three-step staggered inte- 
grator in order to improve accuracy [13], especially when the fluid equations are 
nonlinear and treated implicitly [14]. 

The objective of this paper is the investigation of a broader range of par- 
titioned procedures for the transient solution of coupled aeroelastic problems, 
with particular attention to accuracy and stability issues, subcycling schemes, 
accuracy v.s. speed trade-offs, implementation on heterogeneous computing plat- 
forms, and inter-field as well as intra-field parallel processing. The complete 
three-dimensional aeroelastic problem is difficult to analyze because it mixes lin- 
ear and nonlinear operators, symmetric and unsymmetric matrices, explicit and 
implicit coupling, and can become physically unstable. Therefore, we begin our 
investigation with the design and analysis of partitioned integrators for a sim- 
plified one- dimensional aeroelastic problem that turned out to be a good model 
problem for the more complex aeroelastic systems that we wish to gain some 
intuition about. We focus on implicit time-integration schemes for the struc- 
tural field, because the aeroelastic response of a structure is often dominated 
by low frequency dynamics. However, we consider both implicit/implicit and 
explicit/implicit fluid/structure partitioned procedures, with and without non- 
trivial prediction schemes. We discuss the computational and implementation 
aspects of each procedure and contrast their respective merits and shortcomings. 
Finally, we validate all the conclusions drawn from the investigation of the model 
problem with the simulation of the two-dimensional transient aeroelastic response 
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of a flexible panel in a transonic nonlinear Euler flow. 


2. A one-dimensional aeroelastic model problem with an Euler flow 


2.1. The piston problem: ALE formulation and linearization 

As a model problem, we consider the one- dimensional piston depicted in Fig. 
1. The equilibrium state of this coupled system is defined by a uniform pressure 
po inside and outside the piston chamber, a uniform gas density po, a zero flow 
velocity u 0 = 0, and a chamber length equal to Iq. The gas is assumed to be 
perfect, and the flow isentropic. Hence, the pressure p is function of the density 
p only and obeys: 

dp _ 2 

dp 

where c denotes the sound speed. The cross sectional area of the chamber is 
assumed to be constant and equal to one. 


( 2 ) 



Equilibrium length of the chamber 


x-axis 



Fig. 1. The one-dimensional piston problem 


For this model aeroelastic problem, the one dimensional mass and momentum 
conservation equations for the fluid are: 


dp d'. 

at + di {pu) = 0 

i-( P u) + -?-((pu' +p)) = 0 

The linear dynamic equilibrium of the piston is governed by: 

I 

mq + dq + kq = p(l 0 + q) - Po 


( 3 ) 


( 4 ) 
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where m, d, k , and q denote respectively the piston mass, damping, stiffness and 
displacement. A dot superscript designates a derivative with respect to time. 

The boundary conditions for this coupled fluid/structure problem are given 

u(l 0 + q) = q 


Eqs. (5) above state that the fluid velocity is zero at the fixed wall, and equal to 
the piston speed at the other end of the chamber. 

Clearly, the fluid flow has one moving boundary. Therefore, it is convenient 
to re-write Eqs. (3) with respect to a moving frame characterized by a velocity £ 
that may be different from the fluid velocity u and from zero. Let J = det(dx/d £ ) 
denote the jacobian of the frame transformation x — * f . The ALE form of Eqs. 
(3) goes as follows: 

= 0 

J^ J P U ^ + §^^P U ( U -Z)+P)) = 0 

The above equations can be re-written in vector form as: 

-° < 7 > 



where W and F c are respectively the fluid state and fluid convected flux vectors: 


F c = 


P(u ~0 

(pu(u-i)+p) 


) 


( 8 ) 


The convection matrix associated with the above convected flux vector is: 


m = 


dF c 

dW 


= ( -i i 

y — u 2 + c 2 2u — £ J 


( 9 ) 


We consider the response of the aeroelastic coupled system to small pertur- 
bations around the equilibrium position ( po , «o = 0, po, c<j). First, we note that 


the fluid state vector at equilibrium W( 


•-M 


satisfies: 


j§ (JWa) + i (F ‘ (Wa)) - 0 


( 10 ) 
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and the convection matrix at equilibrium is: 


*«» - (cl o) (“) 

Then, we linearize the convected flux vector around Wo: 

F C (JF) = F C (1F 0 ) + J 0 (0)6W - 8£W 0 (12) 

Finally, from Eqs. (7-12) it follows that the linearized fluid flow equations are: 

~(JSW) + -t(U0)6W - 6(W a ) = 0 (13) 


2.2. Spatial discretization: finite volume formulation and upwinding 


The one- dimensional chamber region is discretized into N grid points and N 
cells (Fig. 2). Integrating Eqs. (13) between Xj_| and x J+ ^ and using a finite 
volume formulation with upwinding leads to: 


Ax, SW + iF^-F^) 

where 


= 0 


Ax,- = 


J j-n — Sj-l 

2 


(14) 


Here, F? +i and F c _^ are the convected numerical fluxes [7] associated with the 
classical +/— flux splitting: 


*7+1 = Jtmswj + J 0 -(0){W, + , - S( i+i {W Zf .t w Hv) 

F'-i = J}(0)SWj-, + J o -(0 ]SWi - W °‘ } 


and Jq(0 ) and J 0 (0) satisfy J o (0) = Jq'(O) + J 0 (0) and are given by: 


4(0) = 5 
4(0) = \ 



(15) 


(16) 
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_cell 

• • • — h — •— H — • « • 

j-1 “ j j+1 


Fig. 2. Discretization of the one-dimensional flow 
Substituting Eqs. (15) into Eqs. (14) gives: 

A Xj SW - J+(0)6Wj. x + ( J 0 + ( 0) - J 0 -(0))^- + J^(0)8W ]+1 




i+ 


(Wo, +W 0j+1 ) 


s tj-, 


(Wbi-t+Wp,) _ 


= 0 


(17) 


2.3. Transpiration 

From Fig. 2 and the second of Eqs. (5), it follows that: 

S( N+i = q (18) 

All other ALE grid velocity perturbations 8£j, j = 1, N axe arbitrary. In order 
to simplify the piston problem from a three- to a two-field coupled problem, we 
assume that these velocity perturbations are small compared to the unperturbed 
sound speed Co- Consequently, Eqs. (17) become: 

Axj 8W - J+WSWj-x + (J 0 + ( 0) - Ji(0))6W j + Jo(0)6W j + 1 =0 (j ± N ) 
Ax n 8W - J+(0)8W n .x + (4(0) - Jq(Q))8W n - qW 0 „ = 0 

(19) 

The quantity qWo N = corresponds to a “transpiration” flux. The reader 

can check that except for the presence of this transpiration flux, Eqs. (19) axe 
identical to the semi-discrete linearized equations governing a one-dimensional 
fluid flow with fixed boundaries. 
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( 20 ) 


2-4- The semi- discrete aeroelastic model problem 
We define the structure state vector as: 

«-(!) 

Using Q, the structural Eqs. (4) can be re-written as: 

Q = ( —Jl i ^ ( p( f o+q)-po ^ (21) 

\ m ms \ m / 

For the linearized piston problem, the forcing term ^ r<'*- | -T)-f° ) ™(.Spn ,) 

= (%u))- Hence ’ E<1 ' (21) 


is a linear function of the fluid state vector SW 
can be re-written as: 

Q = DQ + CSW 
where 


D = (A A) 


( 22 ) 


Also, Eqs. (19) can be re-arranged in matrix form as follows: 

SW = ASW + BQ 

where B is the matrix induced by the transpiration flux qW $ N . 


(23) 


In summary, the semi-discrete coupled system associated with the one- 
dimensional aeroelastic model problem is completely defined by: 


(SW\ 

V Q ) " 

(A B\ 

\C DJ\Q ) 

/<5W\ 

(sw\ 

v Q /(<= o) 

V Q Jo 


(24) 


In the remainder of this paper, we focus on developing, analyzing, and val- 
idating partitioned procedures for solving Eqs. (24). Because the aeroelastic 
response of a structure is often dominated by low frequency dynamics, we con- 
sider only implicit schemes for time-integrating the structural field. However, 
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we consider both explicit and implicit time-integrators for advancing the fluid 
field, as both approaches are popular in computational fluid dynamics. Elegant 
methods for analyzing the stability of partitioned integrators with and without 
subcycling can be found in [15,16]. However, both of these references deal with 
symmetric fields only. We have found that the extension of these analysis meth- 
ods to mixed symmetric/unsymmetric problems such as those described by Eqs. 
(24) is difficult — if not impossible — which has also motivated us to investigate 
first a simplified aeroelastic model problem. 

REMARK 1. Eqs. (24) are also valid for two-dimensional and three-dimensional 
linearized aeroelastic problems. 


3. Mathematical preliminaries 

3.1. Physical v.s. numerical instabilities 

Transient fluid (gas)/structure interaction problems have one particularity: 
they possess a wide variety of self-excited vibrations and instabilities. For exam- 
ple, at speeds of flow somewhat above the critical flutter speed [17], the structural 
system extracts energy from the flow system and a small accidental disturbance 
of the airfoil can serve as a trigger to initiate an oscillation of great violence. 
Physical instabilities can also occur in the linear regime. An example of a linear 
dynamic instability is vibrations due to von Karman vortices [18]. If the frequency 
of the structure loading caused by the vortices is close or equal to the natural 
frequency of the body, then a resonance effect is present and large amplitudes of 
vibrations result. Therefore, when it comes to analyzing the numerical stability 
of a proposed algorithm for time-integrating fluid/structure interaction problems, 
it is essential to consider the case where the coupled system is physically stable 
— that is, when Eqs. (1) or even Eqs. (24) have a solution that does not grow 
indefinitely in time. 

The objectives of this section are to present a mathematical framework for 
the stability analysis of the solution of the semi-discrete Eqs. (24), and to show 
that for the aeroelastic model problem introduced and discretized in Section 2, 
these equations have always a stable solution. Hence, the fluid/structure interac- 
tion model problem presented in this paper is also a good problem for analyzing 
partitioned time-integrators with particular reference to numerical stability. 

REMARK 2. Intuitively, one can expect Eqs. (24) to admit a stable solution for 
the aeroelastic model problem, because the fluid flow is confined inside a closed 
chamber and therefore has a limited amount of energy to exchange with the 
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piston, and the piston is not excited by any other external and t im e-dependent 
force. However, the analysis framework presented in Section 3.2 is interesting 
because it also reveals a numerical property of the coupled model problem that 
turns out to be important for the design of an unconditionally stable partitioned 
procedure for solving Eqs. (24). 

3.2. Analysis framework 

Let X , M, and Sp ( M ) denote respectively a real vector, a real matrix that 
is diagonalizable in the complex space C, and its set of complex eigenvalues. We 
focus our attention on the linear system of ordinary differential equations (ODE): 

X = MX (25) 

First, we introduce two definitions. 

DEFINITION 1. We will say that M is “stable” if and only if: 

a) M is diagonalizable in C 

b) VA € Sp (M), SR (A) < 0 

In b), (A) designates the real part of A. 

DEFINITION 2. We will say that the real symmetric positive definite 
(RSFD) matrix Em is an “energy matrix” for M if and only if EmM is non- 
positive, that is: 

VX, X*E m MX < 0 

Here, the superscript t designates the transpose operation. 

Next, we state and prove four theorems. 

THEOREM 1. An RSPD matrix Em is an energy matrix for M if and only 
if: 

VX solution of X = MX, ^-(^X^mX) < 0 (26) 

dir A 

Proof. From Eq. (25), it follows that j i {\X t E M X) = X*E m MX. Hence, 
■j^{^X 1 'EmX) < 0 if and only if Em is an energy matrix for M. 

THEOREM 2. If M = P 1 HR denotes the diagonalization of a stable matrix 
M, then an energy matrix for M is given by: 

E m = P'P + P*P (27) 
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where P is the complex conjugate matrix of P. 

Proof. Clearly, Em — P* P + P*P is read symmetric. For all real vectors 
X , X l E\fX = 2|PX| 2 is positive and equal to zero only if X = 0 because P is 
non singular. Finally since M = M, we have: 

X'EmMX = X'P'SIPX + X'P'Pp-'toPX 

= 23? (X'P'SIPX) =2j^|(PX) i | 2 3?(f2 t ) < 0 


which completes the proof of THEOREM 2. 

THEOREM 3 (reciprocal of THEOREM 2). Let M be a real matrix that is 
diagonalizable in C. If there exists an energy matrix Em for M, then M is stable. 

Proof. Let X = R + il ^ 0, where i 2 — —1, denote a complex eigenvector 
of M associated with an eigenvalue A. If Em is an energy matrix for M, we have: 

0 > R'EmMR = R*Em*{MX) = R*E M 3t( AX) = R 1 Em[% (A)P - 5 (A)/] 
0 > I 1 E M MI = I* E m Z (MX) = PE M Qt( \X) = I'EmI® (A)J + Of (A)P] 

(28) 

where S (A) designates the imaginary part of A. Adding the two inequalities in 
(28) and exploiting the symmetry of Em leads to 3? (A)[P'Pa*P + PEmI] < 0. 
Since Em is RSPD, it follows that 3 ? (A) < 0, which completes the proof of 
THEOREM 3. 

THEOREM 4. If A and D are two real stable matrices with energy matrices 
Ea and Ed, then: 


E a B + (EdCY = 0 


M 


-a s)' 


is a stable matrix 


Proof. The matrix E m = ^ 
vector satisfying X = MX , we have: 


is RSPD 


. If X = ( ^ is a real 


Q ) 


j^X'EmX) = SW*E a A6W + SW'EaBQ + Q*E D C6W + Q'EdDQ 

= 6W*E a A6W + 8W\ [e a B + (EdC)*] Q + Q'EdDQ 
= SW'EaASW + Q t E D DQ < 0 
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> 


which in view of THEOREM 1 implies that Em is an energy matrix for M. From 
THEOREM 3, it follows that M is stable. 

Theorems 1-3 set the stage to THEOREM 4 which has a nice physical in- 
terpretation. An uncoupled fluid system is physically stable: it does not produce 
energy. If D is non-positive, an uncoupled structural system is also stable. For 
a coupled fluid/structure system, EaB + ( EdC)* = 0 simply expresses that the 
energy extracted from one system is equal to that injected into the other one. 
Hence, THEOREM 4 merely states that a coupled system where the energy is 
exactly conserved is a physically stable system. 

3.3. Physical stability of the model problem , 


Consider again the aeroelastic model problem introduced in Section 2. As- 
suming a constant mesh size Ax, the linearized energy £ fluid °f the discretized 
fluid system can be written as: 


£,uu = Ax 


) 


(29) 


U V« 2 

and its perturbed state vector is SW = (6p\, 6 (pqUi), . . . , Sp^, S(p 0 u^)) . There- 
fore, Ea can be constructed for this system as follows! 


/ 


E a = 


P o 

0 


0 

Ax 

Po 


V o 


Po 

0 


0 \ 


0 

J 

Po ' 


(30) 


Using Eqs. (19), the reader can verify through tedious but elementary calculations 
that Ea is an energy matrix for A induced by the spatial flux splitting. 

For the piston, the state vector is Q — and the energy is £ p ,s* on = 

( kq 2 + mq 2 )/ 2. Hence, for this structural system Ed can be written as: 


= (o m) 

Proving that Ed is an energy matrix for D is straightforward. Using the second 
of Eqs. (22) we have: 

fk 0 \ / 0 • 1 \ _ / 0 k\ 

VO -J W ~ \-k -d) 


E d D = 


12 


which shows that EdD is negative when the damping d is positive, and therefore 
proves that Ed is an energy matrix for D. 


The transpiration term in the last grid cell and the pressure force on the 
piston generate respectively the matrices: 


B = 


(l 

0 

... 0 

° ^ and C = 

(° 

0 ••• 

0 

J2 


v° 

0 

. . . -py 

Ax 

0 ) 


0 ••• 

2sl 

m 

0; 


(32) 


From Eqs. (30-32) and after some algebraic manipulations, it follows that 
EaB-\-(EdC) 1 = 0, which shows that the semi-discrete aeroelastic model problem 
introduced in Section 2 admits a stable solution. Therefore, staggered algorithms 
for time-integrating Eqs. (24) can be analyzed for unconditional stability. 


4. A family of implicit/implicit partitioned procedures 

4-1. Unconditionally stable staggered time-integrators 

Here, we present a family of unconditionally stable implicit/implicit stag- 
gered algorithms for solving the model Eqs. (24) whose “design” is based on the 
following 4-step methodology. 

Step III. Predict the structural field using the value computed at t n = nAt: 

Q p = Q n 

Step 112. Advance the fluid system using the trapezoidal rule: 

* = SW n + A tA6W n+ $ + AtBQP 

where 

,1 SW n + SW n+1 

,SW n +5 = -purr 

2 

Step 113. Advance the structural system using an implicit time-integrator se- 
lected from the so-called generalized trapezoidal family of methods [19], 
and a midpoint value of the previously updated fluid state SW n+ 

Q n+1 = Q n + AtDQ^ + AtC6W n+ * 
where 

Q n+a = (1 -a)Q n +aQ n+1 a € ]0, 1] 
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Step 114. Correct the equations giving 6W n+1 and Q" +1 to enforce uncondi- 
tional stability of the implicit/implicit staggered procedure: 


SW n+1 = <5W n+1 + [6W n+lC ] 
Q n+ 1 = cF^ + [Q n+1 1 


It should be emphasized that the above steps describe the design process of a 
solution methodology and not the computer implementation of a time-integration 
algorithm. In particular, neither the fluid nor the structural fields should be solved 
until the correction terms [<SW n+1 ] and [Q n+1 ] are first specified. 

The trapezoidal rule is unconditionally stable a nd second-order accurate 
when applied to the solution of the uncoupled linearized fluid system (Eqs. (17)). 
For a > 1/2, the generalized trapezoidal family of methods is unconditionally 
stable when applied to time- integrate the uncoupled structural problem. For 
a € ]0, 1], these methods are first-order accurate, except for a = 1/2, in which 
case the corresponding scheme is second-order accurate. The correction terms 
[6W n+1<: ] and [Q n+1 °] should be computed to ensure the unconditioned stability 
of the resulting implicit/implicit staggered solution procedure. 

THEOREM 5. For a > 1/2, [6W n+lC ] = §A t 2 BC8W n+ > and [Q n+1 ‘] = 
(1 - a)At 2 DC6W n+ * , the implicit/implicit staggered time-integrator defined by 
Steps II1-II4 above is unconditionally stable. 

Proof. For [<51¥ n+lC ] = \A t 2 BC6W n+ * and [Q n+lC ] = (l-a)At 2 DC6W n+ $ , 
the proposed implicit /implicit staggered solution edgorithm for solving Eqs. (24) 
becomes: 


(33) 


SW ” +1 = 6W + AiA6W n+ l + AtBQ" + 1a t 2 BCSW n+ i 

£ 

Q n+1 = Q n + AtDQ n+Q + AfCW n+ » + (l - a)At 2 DC8W n+ i 

« e [i l] 


The above partitioned procedure can also be written as: 

Q* = Q n + ^C8W n+ ± 

£ 

Q** = Q n + AtC6W n+ ± 

6vr n +i _ sw n ± . 

— = A8W n+ * + BQ* 

At 


( 34 ) 


Q n+1 - Q* 

At 


= D((l-a)CT +aQ n+1 ) 
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Using the energy matrices E& and Ed, we define the system energy as: 


Sw,Q = lsw t E A SW + ±Q t E D Q (35) 

From Eqs. (33-35) and DEFINITION 2 and after some algebraic manipulations, 
it follows that: 


S 6Wn+ i iQn = S 6W n iQn + At6W n+ *‘ E A A8W n+ > + AtQ*' B'E A 8W n +± 

— < £sw n + 1 ,Q n ^ £sw n ,Q n + A tQ* B t E A 8W n 

A i 2 

£W»+i,q*» = £gW”+\Q” + At8W n +*' C'EdQ" + =f-8W n+ * t C'EdCSW 71 ** 

Cd 

— ^ £sw n + i ,Q •• = £gw n + l ,Q n + A C 1 EdQ* 

£ 6 wn+\Q »+ 1 = £sw*+W + A<((1 - a)Q** + aQ n+1 )‘E D D((l - a)Q** + aQ n+1 ) 

/\/2 

+ (1 -2a)— ((1 “«)<?** +<*Q n+1 ) t D t E D D((l — ot)Q** + aQ n+1 ) 

=S> For a > -, ^V^" + l,Qn + l < £gwn + l } Q*» 

(36) 


which also implies that: 


£gw n + l ,Q -+ 1 < £W« ,Q» + AtQ*' B < E A <5IF n+ ^ + A*5JF n+ i' ‘ C'EdQ* 

=* *W-+i ,Q»+1 < £s W n tQn + A t8W n+ >‘ [E a B + ( E d CY)Q * 

(37) 

Finally, since the aeroelastic model problem satisfies E A B + ( EoCf = 0, it 
follows that: 

£«W"+ 1 ,Q"+i < £gw n ,Q » (38) 

which shows that the numerical energy of the system does not increase in time, 
and therefore the partitioned solution procedure (33) is unconditionally stable. 

THEOREM 6. The implicit/implicit unconditionally stable partitioned procedure 
defined by Eqs. (33) is first-order accurate. 

Proof. Expanding the various terms in Eqs. (33) around the time nAt leads 
to: 

8iv n = A8W n + BQ n + O (At) 

Q" = C8W n + DQ n + O (At) 
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Comparing the above equations with Eqs. (24) completes the proof of this theo- 
rem. Clearly, second-order accuracy would require a more sophisticated predictor 
in Step III. 

It is interesting to note that with the partitioned solution methodology de- 
scribed in Steps II1-II4, we are able to achieve unconditional stability without 
resorting to an augmentation technique [8,20]. Augmentation based schemes are 
often expensive and cumbersome to implement because they require forming and 
factoring the product of the independent field and coupling matrices. The stag- 
gered time-integrator described in Eqs. (33) requires only one additional sparse 
matrix-matrix product to form BC. 

REMARK 3. Under some mild assumptions, the condition EaB + {EqC) 1 = 
0 can be shown to hold for two-dimensional and three-dimensional linearized 
aeroelastic problems. In that case, the staggered time-integrator described in 
Eqs. (33) is also unconditionally stable for these two-dimensional and three- 
dimensional problems. 

4-2. Subcycling 

The fluid and structure fields have often different time scales. For problems 
in aeroelasticity, the fluid flow usually requires a smaller temporal resolution 
than the structural vibration. Therefore, if the unconditionally stable staggered 
algorithm (33) is used to solve a coupled fluid/structure problem, the coupling 
time-step A t c will be typically dictated by the time-step A tp that guarantees a 
certain accuracy in the flow solution, rather than the time-step A t$ > A tp that 
meets the accuracy requirements of the structural field. 

Using the same time-step A t c in both fluid and structure computational 
kernels presents only minor implementational advantages. On the other hand, 
subcycling the fluid computations with a factor ns/F — Ats/ A <f can offer sub- 
stantial computational advantages, including: 

• savings in the overall simulation CPU time, because in that case the struc- 
tural field will be advanced fewer times. 

• savings in I/O transfers and/or communication costs when computing on a 
heterogeneous platform, because in that case the fluid and structure kernels 
will exchange information fewer times. 

However, the computational advantages highlighted above are effective only 
if subcycling does not restrict the stability region of the staggered algorithm 
to values of the coupling time-step A t that are small enough to offset these 
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advantages. For example, consider the following fluid-subcycled version of the 
unconditionally stable staggered time-integrator given in Eqs. (33): 

w +1<0) = SW n 

{ 

For k = 0, ..., ns/F — 1 

(5 ^+i (i+1) _ sW n+l(k) + (AtA + ^-BC)SW n+l(k+ ^ ) + AtBQ n (39) 

} 

^Tl + l _ ^"+1 ( " 5/F) 

Q n+1 = Q n + A tDQ n+a + A tC8W n +± + (1 - a)At 2 DC6W n+ * 

This algorithm implements the simplest possible subcycling scheme and is of- 
ten used in many applications. Unfortunately, the reader can easily check that 
the above fluid-subcycled partitioned procedure (39) is no longer unconditionally 
stable. Next, we present an improved subcycling approach that preserves the 
unconditional stability of the partitioned procedure (33). 

THEOREM 7. For a > 1/2, the following fluid-subcycled version of the 
staggered time-integrator given in Eqs. (33) is unconditionally stable: 

^n+l (0) _ 8W n 

= Q n 

{ 

For k = 0, ..., ns/F ~ 1 

< 5 ^n+i (l+1) = sw n+l(k) + ( AtA + BC)SW n+l(k+ + A tBX {k) 

X (fc+1) = X {k) + A<CW n+l( ‘ + * > 

} 

6W n+1 = 6JU n+l( " s/F) 

Q n+ 1 = + n s/F A*D((l -a)X (ns "> +aQ n+1 ) 

* e i u 

(40) 

Proof. The proof of this theorem is siipilar to that of THEOREM 5 and uses 
the system energy defined in Eq. (35). Using Eqs. (35,40) and DEFINITION 2 
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one obtains: 1 

£*w"+i (0) ,x<°) = €sw«,Q” .... ,- ; , ur . 

e »w - »<>+»,*<« ^ + Ai«W'" +1 “ +1> ‘.E.,B(.t<*> + ^CW> +,< * H> ) 

^ ^ e tw , w . ) , xw +&t6W"+ ll ‘+ i> ‘E A BX< k +» 

£ w«<‘+'>, *<>■»> = *W« + .“+‘>,j<m + AfW” +I< * +1 ’V,MX<‘> + ^CW +1 “ +t 

=*• £ w . t .<‘+» ,*»♦.> = £,„,„+.(*«> ,x<‘) + A«W"* +1 “ + *’‘ C‘E D X^ 

£sW n + l ,Q n+l = ^ 6VV „+i (n s/i : ’) (Xn+1 ^s/f) 

+ n s/F At((l - a)X^ ns / F) + aQ n+1 YE D D((l - a)X ins / F) + <*Q n+1 ) 

+ (1 - 2c*)^((l - a)X (ns /^ + aQ n+1 )* D 1 E d D((1 - a)X (ns /^ + aQ n+1 ) 

— '< For o > — , ^giv n +* ,Q n+l — ^(y( # s/F) ^( b s/f) 

(41) 

The above inequalities also imply that: 

k=n s/F -l t 

^»r.+i,Q-+i < Ssw^qn + At Y, SW n+l(k+i) [E A B + (E D cy]X( k+ » 

k-0 

(42) 

Finally, since the aeroelastic model problem satisfies EaB -f- (EdC)* = 0, it 
follows that: 

£ew n + 1 ,Q n + l < £sw n ,Q n (43) 

which shows that the numerical energy of the system does not increase in time, 
and therefore the partitioned and fluid-subcycled solution procedure (40) is un- 
conditionally stable. 

THEOREM 8. The implicit/implicit unconditionally stable partitioned and fluid- 
subcycled procedure defined by Eqs. (40) is first-order accurate. 

Proof. The proof of this theorem is similar to that of THEOREM 6. 

The subcycling approach advocated in Eqs. (40) preserves the computational 
advantages of subcycling. At each stage, the evaluation by the fluid solver of 
the correction term does not require neither advancing the structural 

state vector, nor exchanging information with the structural solver. Moreover, 


updating and BX^ requires only two sparse matrix-vector products and 

therefore is relatively inexpensive. 

The interpretation of the role of the correction term goes as fol- 

lows. In order to solve Q = DQ + C6W between t n and t n +\, the struc- 
ture kernel must receive from the fluid module the best possible approxima- 
tion of the coupling quantity ff” +1 CSWdt. In the fluid-subcycled partitioned 

procedure (39), this integral is approximated by A tC8W n +%, which guaran- 
tees a certain accuracy but does not warrant unconditional stability. On the 
other hand, the strategy consisting in approximating CSWdt via updat- 

ing = X + A tCSW n+1< ' and replacing in the procedure (39) the 

“frozen” A tBQ n by the updated quantity A tBX^ not only provides a better 
coupling accuracy, but also preserves the unconditional stability of the original 
non subcycled partitioned procedure (33). 

The implications of the above results and discussion on the staggered and 
subcycled solution of more complex aeroelastic problems can be formulated as 
follows. When the fluid field is subcycled and updated ahead of the structural 
field, then: 

• the motion of the moving fluid boundary induced by the structural deforma- 
tion should not be completely absorbed during the first fluid sub cycle and 
“frozen” during the remaining ones. Rather, this induced motion should be 
distributed among all subcycle stages via a careful interpolation scheme. 

• after all fluid subcycles axe completed, the mean value rather than the final 
value of the pressure field must be transmitted to the structure. 

4-3. Examples 


Here, we illustrate the numerical properties of the family of implicit/implicit 
staggered procedures presented in Sections 4. 1-4.2 with the solution of the aeroe- 
lastic model problem (24). We consider the case where or = 1/2 and the structure 
is undamped (d = 0). First, we introduce the non-dimensional variables: 



( 44 ) 
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and rewrite Eqs. (19) and (4) in non-dimensional form as follows: 


A xj SW'j - X (0)<5 VTy-! + (7; (0) - J 0 (0))t> W, + J 0 (0 )SW j+1 =0 (j ^ N) 

Ax N ^V-7j(0)6F^_ 1 + (J^(0)-7r(0))^W : N -( 9 0 ^ = o 


9 +w‘g = 
where 

t 



u ) 2 


m — 


- !) 
m 


d_ 

di 

1(5 5)’ ^ (0) - K'i 1 -i) 

i± 1 

c§ m 
m 

Po^o 


(45) 


Next, we discretize the piston chamber into 21 grid points and N = 20 finite 
volume cells, and set the non-dimensional parameters to u> 2 = 3.03 x 10 -4 and 
m = 30.77. We consider the following initial conditions: 


Sp{t = 0) = Spu(t = 0) = 0 

q(t = 0) = 0 q(t = 0) = 1 


and solve the coupled equations (45) using eight time-steps varying between At = 
lxCFL and At = 128xCFL, where CFL = lo/(Nc 0 ) is computed with respect to 
the uncoupled fluid problem. The obtained non-dimensional piston displacement 
q/lo and fluid pressure in the cell in contact with the piston (p/ pocfyio are depicted 
in Fig. 3-4 for the case without subcycling. 

Clearly, the results reported in Fig. 3-4 highlight the unconditional stability 
of the family of implicit/implicit staggered procedures presented in Sections 4.1- 
4.2. Stable responses are observed for all time-steps, and accurate results are 
obtained for both the piston displacement and fluid pressure for time-steps as 
large as At = 32 x CFL. 

The previous computations are repeated using the fluid-sub cycled staggered 
time-integrator (40), 2xCFL < Atp < 4xCFL, and several subcycling factors 
1 < ng/F 32. The corresponding results (Fig. 5-8) confirm numerically the 
unconditional stability proved mathematically in THEOREM 7. Note however 
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that large subcycling factors ns/F introduce a spurious phase shift in the initial 
stages of the coupled computations that can ruin the accuracy of the response 
history. Both amplitude and phase errors can also be observed in that case. This 
suggests that an adaptive time-stepping strategy is needed in order to resolve 
better the initial response of the coupled system. 









Fig. 4. Computed non-dimensional fluid pressure (without subcycling) 












II Scheme with subcycling - di=2xCFL n Scheme with subcycling - di=2xCFL 



0 82 165 248 330 0 82 165 248 330 

— Ns/F=4 — Ns/T=8 

Fig. 6. Computed non-dimensional fluid pressure (At = 2xCFL and with subcycling) 


25 











II Scheme without subcycling - dt=4xCFL 


II Scheme with subcyding - dt=4xCFL 



Fig. 8. Computed non-dimensional fluid pressure (A t = 4xCFL and with subcycling) 






5. A family of explicit/implicit fluid/structure partitioned procedures 

Next, we consider the case where an explicit scheme is desired for advanc- 
ing the fluid field. Consequently, we cannot design an unconditionally stable 
staggered algorithm for solving the coupled fluid/structure problem, because in 
that case it would be only conditionally consistent [21]. Rather, we focus on 
developing a family of partitioned procedures whose stability limit is governed 
by the critical time-step of the fluid solver. In other words, we wish to design 
a staggered solution algorithm for the coupled problem whose stability limit is 
not worse than that of the underlying fluid explicit time-integrator. This is not 
necessarily a trivial task because coupling effects can restrict the stability limits 
of the independent field time-integrators. 

5.1. A predictor- corrector approach 

Here, we present a family of explicit/implicit staggered algorithms for solving 
the fluid/structure Eqs. (24) that is based on similar ideas to those presented in 
Section 4: 

Step Ell. 


Step EI2. 


Step EI3. 


Step EI4. 

6W n+1 = 5W P + [6W n+lC ] 
Q n+ 1 = + [<9 n+1 1 


Predict the structural field using the value computed at t„ = nAt: 

Q p = Q n 

Predict the fluid system using the forward Euler explicit scheme: 

5W P = 5W n + AtA6W n + AtBQ p 

Improve the structural field using an implicit generalized trapezoidal 
method and the predicted fluid state 6W P : 

QZ+i = Q n + AtDQ*+° + A tCSW p a 6 ]0, 1] 

Correct the expressions yielding the fluid and structural fields to en- 
hance the stability of the explicit /implicit staggered procedure: 


When applied to the solution of the uncoupled linearized fluid system (Eqs. 
(17)), the forward Euler algorithm is first-order accurate and stable for CFL < 1. 
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For a > 1/2, the generalized trapezoidal family of methods is unconditionally 
stable when applied to time-integrate the uncoupled structural problem. For 
a E ]0, 1], these methods are first-order accurate, except for o = 1/2, in which 
case the corresponding scheme is second-order accurate. The correction terms 
[£W n+lC ] and [Q n+1<: ] should be computed to enhance the stability of the resulting 
explicit/implicit staggered solution procedure. 

THEOREM 9. For a > 1/2, [8W n +' c ] = aAtB(Q n+1 - Q n ) and [Q n+1 '] = 
— AiC(AtBQ n + 6W n+1 )/2, the stability of the explicit/implicit staggered time- 
integrator defined by Steps EI1-EI4 above is governed by the stability of the 
explicit time-integrator of the uncoupled fluid problem. 

Proof. For [6fF n+lC ] = aAtB(Q n+1 - Q n ) and[Q rt+1 ‘] = -AtC(AtBQ n + 
SW n+l )/2, the proposed explicit /implicit staggered solution algorithm for solv- 
ing Eqs. (24) becomes: 


Q n+ 1 = Q n + AtDQ n+a + AtC(8W n + AtASW n + —BQ n+Q ) 

2 

SW" +t = SW" + AtASW” + AtBQ n+< ‘ 


(46) 


Let E n be defined as follows: 


E n = | Q nt E D Q n 

If the second of Eqs. (46) is re-written as: 

Q n+l = Q n + AtDQ n+a + AtZ 
where 

Z = C(6W n +AtA8W n + ^-BQ n ) 
it follows that for a > ^ we have: 

£«+i < E n +AtQ n+at E D Z 

Using the above inequality and calculations similar to those in (36), one can 
show that if the energy of the uncoupled fluid problem decreases in time, then 
the energy of the coupled system also decreases in time. Hence, the stability of 
the explicit /implicit staggered procedure defined by Eqs. (46) is governed by the 
CFL condition (CFL < 1) of the explicit Scheme applied to the uncoupled fluid 
problem. 
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THEOREM 10. The explicit /implicit partitioned procedure defined by Eqs. 
(46) is first-order accurate. 

Proof. The proof of this theorem is similar to that of THEOREM 6. 

REMARK 2. Note that while the family of implicit/implicit partitioned pro- 
cedures (33) and its fluid-subcycled version (40) require updating the fluid flow 
ahead of the structural field, the family of explicit/implicit staggered algorithms 
(46) require updating the structural system first. 

5.2. Examples 

Here, we illustrate the numerical properties of the family of explicit/implicit 
staggered procedures presented in Section 5.1 with the solution of the non- 
dimensional coupled equations (45). We consider the case a = 1/2, and use 
the same non-dimensional parameters and finite volume mesh as in Section 4.3. 

First, we solve the coupled problem for 0.5xCFL < A t < l.OxCFL, where 
the CFL is with respect to the uncoupled fluid problem. The obtained non- 
dimensional piston displacement q/lo and fluid pressure in the cell in contact 
with the piston (p/ p 0 cjj)20 are reported in Fig. 9-10 for the case without subcy- 
cling. Clearly, these results demonstrate numerical stability for A t < 1.0 x CFL. 
However, they also show that at time-steps close to the uncoupled fluid CFL con- 
dition, the errors introduced in the initial stages of the computations propagate 
throughout the entire history response of the fluid pressure, but do not affect 
significantly the evaluation of the structural displacement. 

Next, we repeat the previous explicit/implicit computations and subcycle 
the fluid system. We use a subcycling scheme similar to that introduced in Eqs. 
(40) in order to maximize the coupled stability time-step. In that case, the nu- 
merical results reported in Fig. 11 for the non-dimensional piston displacement 
q/k indicate that there exists a maximum subcycling factor beyond which the 
explicit /implicit time-integrator (46) with subcycling becomes numerically un- 
stable. 


30 









El Scheme without subcycling • dt=0.90*CFL 


El Scheme with subcycling - dt=0.90xCFL 



Fig. 11. Computed non-dimensional piston displacement (A t = 0.9xCFL and with subcycling) 







6. Parallel staggering strategies, error analysis, and CPU distribution 


The family of partitioned procedures presented in Sections 4 and 5 are in- 
herently sequential: in the implicit/implicit case, the fluid state vector must be 
updated before the structural system can be advanced, and in the explicit/implicit 
case, the new structural displacements and velocities must be computed before 
the new fluid state vector can be evaluated. With the advent of parallel pro- 
cessors and distributed computing platforms, it becomes interesting not only to 
parallelize each field computations, but also to design staggered time-integration 
algorithms that promote inter-field parallelism — that is, that allow advancing 
simultaneously the fluid and structural systems. In this Section, we present such 
partitioning procedui :s and discuss their accuracy v.s. speed trade-offs. We con- 
sider only explicit/implicit algorithms. * More specifically, we focus on the case 
where the fluid is advanced using the first-order accurate forward Euler scheme 
and the structure is advanced using the second-order accurate midpoint rule, be- 
cause these algorithms are already available in our large-scale simulation parallel 
software. We use the family of time-integrators (46) as reference, and therefore 
begin with their error analysis. 


We introduce the following nomenclature: 


Of 

Os 

T f 

T s 

CPU F 

CPUs 


: number of floating-point operations in one fluid time-step 
: number of floating-point operations in one structural time-step 
: fluid-to-structure single pass transfer time 
: structure- to-fluid single pass transfer time 
: CPU resource allocated to the fluid kernel 
: CPU resource allocated to the structure kernel 


For every partitioned procedure, we give the resource distribution between the 
fluid and structure kernels for simulations on heterogeneous platforms. 


6.1. Algorithm ALGO: the basic staggered scheme 

In the sequel, we refer to the explicit /implicit time-integrator (46) as the 
basic staggered algorithm ALGO. Let Z be defined as follows: 

2 = (*w) 
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Considering the time-interval [t n = nAt, t n+ „ s/F = (n + n S / F )At], it follows 
from Eqs. (46-47) that: 


Z »+ns,r = [/ + ns/F&t (^ 


C\ 2 At 2 { D*+BC 
A ) +n s/F 2 \BD + AB 


+ ^) + O (Af 3 )]Z 


DC + C^4\ 
BC + A 2 ) 


(48) 


If ERsw and ERq denote respectively the errors in the fluid and structural 
responses after ns/ F time-steps, it follows from (48) that: 


ERfk G0 = -n s/F —A6W»+0(At 3 ) 
A / 2 

ER^ Lao = n s/F — CAQ"+Q(At 3 ) 


(49) 


which shows that ALGO is first-order accurate. Hence, the accuracy of the struc- 
tural computation is first-order even though the midpoint rule applied to the 
uncoupled structural problem is second-order accurate. 

The basic steps of ALGO axe graphically depicted in Fig. 12. The CPU time 
needed to advanced the coupled solution n s/ F At is equal to: 

T coup ud = n s/F (T F + T S + + CPUs ^ ( 50 ) 

If the total amount of CPU resources CPU = CPU F + CPUs is assumed to be 
fixed, T COU pi e d is minimum for: 


CPU F _ jo7 
CPUs ~ V Os 


which gives: 


nnALGO 
1 coupled 


n s/F{T F + Ts + 


0 F -f- Os 4" 2y/U^O 

CPU 


(51) 
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Fig. 12. ALGO: the basic staggered algorithm 

6.2. Algorithm, ALGl: subcycling the fluid system 

A fluid-subcycled version of a slightly modified ALGO where the subcycling 
scheme follows the guidelines of THEOREM 7 is referred to as ALGl in the sequel 
and is given by: 

£ W n+l (0) _ 6W n 

X (0) = Q n 

{ 

For k = 0, tis/f “ 1 

^ n+l( ‘ +l) = 6W n+l(k) + AtA6W n+lik) + AtBQ n , ro . 

X (fc+D = X {k) + AtC6W n+l{k+i) 

} 

6W n+1 = SW n+lins ' F) 

Qn+l = Qn+n s/F _ X (n s/F ) + ns/f , Af £)( ^ - 9^12111 ) 


Essentially, the fluid system is subcycled during ris/F time-steps At, and the 
structural field is advanced in one shot a large time-step equal to ns/F^t using 
the average fluid pressure between t n and t n+ „ s/F . 
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Expanding the various terras in Eqs. (52) around t„ = nAt leads to: 


.Isa 


ERffi' = -"l/pyW" + 0 (A t\ ns/pAt*) 

A i3 

ER% LG1 = n 3 s/F =—(D 3 Q n + 2CSW - 2CBQ n ) + 0 (At 4 , n 2 s/F At 3 ) 

(53) 

which shows that ALGl is also first-order accurate. However, from (49) and (53), 
it follows that subcycling amplifies the fluid errors by the subcycling factor ns/F- 
Measuring the effect of subcycling on the structural errors of ALGO is less trivial: 
in order to keep its computer implementation simple, we have designed ALGl 
as a fluid- sub cycled version of a “slightly modified” rather than the “original” 
ALGO. Consequently, the structural errors grow as O (ns/F A< 2 ) in ALGO, and 
as O (n 3 S / F At 3 ) in ALGl. 

The basic steps of ALGl are graphically depicted in Fig. 13. The CPU time 
needed to advanced the coupled solution ris/F^t is equal to: 


Tcoupled = Tp +Ts + 


n s/fOf 

cptTf 


+ 


Os 


CPUs 


(54) 


For a fixed total amount of CPU resources CPU = CPUf + CPUs, Pcoupled is 
minimum for: 

CPU F /^s/fOf 
CPUs ~ V °s 


which gives: 


(55) 


The comparison of % coupled highlights the computational advantages 

of subcycling. 


rnALGl rp , m , n S/F 0 F + Q S + 2 V^S/F^F^S 

coupled F * * S i C 1 T*TT 
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Fig. 13. The fluid-subcycled ALG1 staggered algorithm 
6.3. Algorithm ALG2: improving the accuracy of ALGl 

In order to improve the accuracy of the fluid solution in ALGl, we introduce 
a computational phase shift between the fluid and structure kernels equal to 
ns/FAt/2. Assuming that 8W n and Q n+ ? are available, the improved subcycled 
explicit/implicit algorithm ALG2 computes 8W n+l and Q n+ * as follows: 


SW n + lW = 8W n 

X (0) = Q n+ l = Q n+1 ¥~ 

{ 

For k = 0, ..., n S /F ~ 1 

<W n+ i (fc+1) = 8W n+l ' k) + AtA8W n+l(h) + AtBQ 71 * 21 ^ 

X (k+1) _ x ik) + AtC8W n+l(k+ ^ 

} 

SW "* 1 = £W rw+1< "® /,) 

Q n+ f = Q n+ 2+ n s/F _ x (ns ' F) +ns /F AtD( - n+3 + , Q n+2+ns/f ) 
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Algorithm ALG2 has the same computational and I/O transfer requirements as 
ALG1. However, its error analysis leads to: 


(57) 

A direct comparison of (49), (53) and (57) shows that ALG2 offers the computa- 
tional advantages of ALGl, and the higher accuracy of ALGO. 

6. 4- Algorithm ALG3: introducing inter-field ■parallelism 



ALGO, ALGl and ALG2 are inherently sequential. In all three algorithms, 
the fluid system must be updated before the structural system can be advanced. 
The following explicit /implicit fluid-subcycled time-integrator ALG3 introduces 
inter-field parallelism in the solution of Eqs. (24): 



rn+l (0) 


= 6W T 


For k = 0, ..., rt s /F — 1 

SW n+ l(fc+1) = <W n+l(i) + A tA6W n+l(k) + AtBQ n 

} 

SW n+1 = <5W n+1<ns/ir) 

qu+1 _ Qn+n s , F _ Qn + + n SfF AtD{ ^- - ^ ) 


Clearly, the fluid and structure kernels can rim in parallel during the time- 
interval [t n , t„+„ s/F ]. Inter-field communication or I/O transfer is needed only 
at the beginning of a time- interval. Expanding the various terms in Eqs. (58) 
around t n = nAt leads to: 

ERtiP’ = O (n| /F Ai 2 ) I 

ER$ La3 = O (n| /F A i 1 ) (59) 

which demonstrates that ALG3 is first-order accurate. However, the above error 
analysis also shows that parallelism in ALG3 is achieved at the expense of ampli- 
fied errors in both the fluid (a factor equal to ng/F with respect to ALG2) and 
structural (a factor equal to l/(n 5 /pA<) with respect to ALG2) systems. 
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The basic steps of ALG3 are graphically depicted in Fig. 14. The parallel 
CPU time needed to advanced the coupled solution ns/pAt is equal to: 


TcoupUd - TF+Ts + max( £pjj F > q pu s ) 


(60) 


For a fixed total amount of CPU resources — for example, a fixed number of 
processors in a parallel machine — T coup i e( i is minimum for: 

CPUf _ tis/fOf 
CPUj “ Os 

Hence, the parallel CPU time associated with ALG3 is: 


T ALGZ rp I m I 

coupled — ±F + J-S + 


ks/fOf + O s 
CPU 


(61) 


which demonstrates the computational advantages of this parallel scheme. 



(^^ructure^) 


n+1 


Fig. 14. The basic parallel subcycled ALG3 algorithm 
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6.5. Algorithm ALG4: improving the accuracy of ALG3 

In order to improve the accuracy of the basic parallel time-integrator ALG3, 
we propose to exchange information between the fluid and structure kernels at 
half-step in the following specific manner (ALG4): 


< 5^n+i (0) = sW n 

{ 

For k = o, - 1 

^n+i(*+ l ) _ ^n+i(‘) + At A6W n+lW + A tBQ n 

} 

8W n+ * = SW n+li ~^~ ) 

O n + O n+1 

Q n+ 1 = Q n + ns /F AtC5W n + n s/F AtD (¥— ^ ) 

} 

For k = n s/F ~ 1 

8 W n+l(k+x) = £W n+l( * ) + AtA8W n+lik) + AtBQ^ 

} 

6W n+ 1 = W +1 ‘” s/f) 

Q n+1 = Q n +n s/F AtC5W n+ ± +n s/F AtD(Q— y^) 


(62) 


The above algorithm ALG4 is illustrated in Fig. 15. The first-half of the com- 
putations is identical to that of ALG3, except that the fluid system is subcycled 
only up to t 


n+ 


*S/F 

2 


while the structure is advanced in one shot up to t n +„ s/f ,. 

At t >» S /f , the fluid and structure kernels exchange pressure, displacement and 
2 

velocity information. In the second-half of the computations, the fluid system is 
subcycled from t , "s/f to t n+ns/F using the new structural information, and the 
structural behavior is re-computed in parallel using the newly received pr essure 
distribution. Note that the first evaluation of the structural state vector Q n + l 
can be interpreted as a predictor. 
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Fig. 15. The improved parallel subcycled ALG4 algorithm 
An error analysis of ALG4 reveals that: 


ER(w C * = O (n s/P A i 2 ) 
ER* lc ' = 0(nJ /f .Af 3 ) 


which shows that this parallel algorithm has the same accuracy as the improved 
ALG2. 

The parallel CPU time needed to advanced the coupled solution ris/F^t 
using ALG4 is equal to: 


2111-0 O 

Tcoupted = 2 (Tp + Ts + max{ j CPUs )) ( 64 ) 

For a fixed total amount of CPU resources, this parallel time is minimum for: 

CPUf ks/fOf 

CPUs = 20s 


Hence, the parallel CPU time corresponding to ALG4 is: 


T* L JZL = 2 (7> + Ts) + 


^s/fCf + 20 3 

CPU 


(65) 
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In summary, the accuracy of the basic parallel algorithm ALG3 is improved at the 
expense of an additional communication step or I/O transfer during each coupled 
cycle. 

6.6. Applications 

The advantages and limitations of ALGO — ALG4 are summarized in Fig. 16 
which contrasts the various computed solutions q/lo of the non-dimensional Eqs. 
(45), using the same non-dimensional parameters and finite volume mesh as in 
Section 4.3, and A = 0.9xCFL. 

ALGO- dp4.9Q.CH. ALGl • dWMfcCTl 




ALG3 • dt«O.90»CFL ALW-cWMQiCFI. 




Fig. 16. Computed non-dimensional piston displacement (ALGO — ALG4) 
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If the ALGO solution is used as reference, the reader can observe that ALGl is 
at the stability limit when tis/f — 5. The accuracy of ALG3 is comparable to 
that of ALGl, but its stability is restricted to n$/F = 2. On the other hand, the 
parallel algorithm ALG4 is shown to have the same accuracy as ALGO even for 
n s/F — 20, which highlights the merits of this improved parallel algorithm. 


7. Aeroelastic response of a flexible panel in transonic nonlinear regime 

Based on the insight gained from the analysis and solution of the coupled 
piston problem, we have extended the algorithms presented in this paper to the 
solution of the three-field coupled formulations summarized in Eqs. (1), and com- 
plex aeroelastic problems. The generalization of ALGO — ALG4, their implemen- 
tation on heterogeneous parallel processors, and the analysis of their performance 
results are discussed in details in a companion paper [22]. Here, we focus on vali- 
dating qualitatively the conclusions drawn from the mathematical and numerical 
investigations of the model problem with the simulation of the two-dimensional 
transient aeroelastic response of a flexible panel in transonic nonlinear regime. 

For a two-dimensional simulation, the panel is represented by its cross section 
that is assumed to have a unit length L = 1, a uniform thickness h = 10 -2 x L, 
and to be clamped at both ends. This rectangular cross section is discretized 
into 300 x 3 plane strain 4-node elements with perfect aspect ratio to avoid mesh 
locking. This fine discretization — which generates 1204 nodes — is not needed 
for accuracy; we have designed this structural mesh only because we were also 
interested in assessing some computational and I/O performance issues. The 
two-dimensional flow domain around the panel is discretized into 2880 triangles 
and 1504 vertices. The free stream Mach number is set to M ^ = 0.8, and a 
slip condition is imposed at the fluid/structure boundary. Because the fluid and 
structural meshes are not compatible at their interfaces (Fig. 17), the “Matcher” 
software [23] is used to transfer the pressure load to the structure, and to transmit 
the structural deformations at the surfaces of the panel to the fluid. 

Initially, a steady-state flow is computed around the panel at = 0.8 (Fig. 
18). Next, this flow is perturbed via an initial displacement of the panel that is 
proportional to its second fundamental mode (Fig. 19), and the subsequent panel 
motion and flow evolution are computed using the ALGO — ALG4 explicit/implicit 
fluid/structure time-integrators (Fig. 20-21). 
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Fig. 17. A partial view of the structure and fluid discretizations 


More specifically, the dynamic equations of equilibrium of the structure are 
solved via the parallel implicit transient FETI method [24], with the improve- 
ments proposed in [25] for the efficient iterative solution of systems with repeated 
right hand sides. The Euler flow equations are solved with a parallel algorithm 


that combines a second-order accurate Monotonic TJpwinding Scheme for Con- 
servation Laws for spatial approximation, and a second-order low-storage explicit 
Runge-Kutta scheme for time- integration [26]. All computations are carried out 
on an iPSC-860 parallel processor. Four processors are allocated to the fluid 
code, and two processors to the structural program. The fluid and structural 
computations are implemented in a heterogeneous manner using the intercube 
communication procedures described in [27]. 

/ } . 



\ 

X \ 

Fig. 18. Pressure isovalues of the steady-state flow solution 


Fig. 19. Initial perturbation of the panel displacement field 
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Fig. 20. Pressure isovalues at t = 8.75 X 10 5 



\ \ 

\ 

v \ 


Fig. 21. Pressure isovalues at t = 1.25 X 10 4 
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For the uncoupled fluid problem, the CFL stability time-step is At = lxCFL 
= 6.25 x 10 -6 . The lift solutions computed by ALGO (At c = A t F = A ts = 6.25 x 
lO" 6 ), ALG1 (A t F = 6.25 x 10" 6 , At s = 5 x 10" 5 , n s/F = 8), ALG3 (At F = 
6.25 x 10 -6 , A ts — 5 x 10 — 5 , ns/ F — 8), and ALG4 (A t F — 6.25 x 10 Ats — 
5 x 10 -5 , n S / F = 8) are depicted in Fig. 22. Clearly, as predicted by the theory 
presented in this paper, all proposed explicit/implicit time-integrators are shown 
to be numerically stable at At F = lx CFL and ns/ F = 8, and ALG4 is reported 
to improve the accuracy of ALG3. Note also that the parallel algorithm ALG4 
is shown to achieve in practice a better accuracy than the sequential algorithm 
ALG1 for the same time-steps A t F = lx CFL and Ats = 8x CFL. 



Fig. 22. Computed lifts (ALGO, ALGl, ALG3, ALG4) 

The superiority of ALG4 over ALG3 is also illustrated in Fig. 23 which 
shows that the lift solution computed by ALG4 with Ats = 16 x CFL is less 
oscillating than that computed by ALG3 with a smaller Ats = 8x CFL. While 
both computed lift solutions may not be accurate enough for structural analysis 
purposes, ALG3 and ALG4 quickly (faster than ALGO) and correctly reproduce 
the overall aeroelastic behavior of the system — for example, they show that 
flutter is not occuring — which is what a designer is mostly interested in verifying 
initially. 
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Fig. 23. Computed lifts (ALGO, ALG3, ALG4) 

8. Closure 

In this paper, we have presented several partitioned procedures for time- 
integrating the transient coupled aeroelastic problem, and have discussed their 
merits in terms of accuracy, stability, heterogeneous computing, I/O transfers, 
subcycling, and parallel processing. All theoretical results have been derived for a 
one-dimensional piston model problem with a compressible flow, because the com- 
plete three-dimensional aeroelastic problem is difficult to analyze mathematically. 
However, the insight gained from the analysis of the coupled piston problem and 
the conclusions drawn from its numerical investigation have been confirmed with 
the numerical simulation of the two-dimensional transient aeroelastic response 
of a flexible panel in a transonic nonlinear Euler flow regime. In particular, we 
hope that with the unconditionally stable implicit-implicit staggered procedure 
and the parallel coupling strategy with superior accuracy properties presented in 
this paper, large-scale transient aeroelastic computations will be finally feasible. 



49 




Acknowledgments 

The second author acknowledges partial support by NASA Lewis Research Cen- 
ter under Grant NAG3-1273 and partial support by the Nat ional Science Foun- 
dation under Grant ASC-9217394. The authors would like to thank P.-S. Chen, 
L. A. Crivelli, S. Lanteri, and N. Maman for the parallel implementations on the 
iPSC-860. — 

References 

[1] J. Donea, An arbitrary Lagrangian-Eulerian finite element method for tran- 
sient fluid-structure interactions, Comput. Meths. Appl. Mech. Engrg. 33 
(1982) 689-723rr' 

[2] T. J. R. Hughes, W. K. Liu and T. K. Zimmermann, Lagrangian-Eulerian 
finite element formulation for incompressible viscous flows, U.S. -Japan Seminar 
on Interdisciplinary Finite Element Analysis, Cornell Univ., Ithaca, NY, Aug. 
7-11 (1978). 

[3] T. Belytschko and J. M. Kennedy, Computer models for subassembly simu- 
lation, Nucl. Eng. Design 49 (1978) 17-38. 

[4] 0. A. Kandil and H. A. Chuang, Unsteady vortex-dominated flows around 
maneuvering wings over a wide range of mach numbers, AIAA Paper No. 88- 
0317, AIAA 26th Aerospace Sciences Meeting, Reno, Nevada, January 11-14, 
1988. 

[5] C. Farhat and T. Y. Lin, Transient aeroelastic computations using multi- 
ple moving frames of reference, AIAA Paper No. 90-3053, AIAA 8th Applied 
Aerodynamics Conference, Portland, Oregon, August 20-22, 1990.. 

[6] J. T. Batina, Unsteady Euler airfoil solutions using unstructured dynamic 
meshes, AIAA Paper No. 89-0115, AIAA 27th Aerospace Sciences Meeting, 
Reno, Nevada, January 9-12, 1989. 

[7] M. Lesoinne and C. Farhat, Stability analysis of dynamic meshes for transient 
aeroelastic computations, AIAA Paper No. 93-3325, 11th AIAA Computational 
Fluid Dynamics Conference, Orlando, Florida, July 6-9, 1993. 

[8] K. C. Park and C. A. Felippa, Partitioned analysis of coupled systems, in: 
Computational Methods for Transient Analysis, T. Belytschko and T. J. R. 
Hughes, Eds., North-Holland Pub. Co. (1983) 157-219. 

[9] C. J. Borland and D. P. Rizzetta, Nonlinear transonic flutter analysis, AIAA 
Journal (1982) 1606-1615. 


50 


[10] V. Shankar and H. Ide, Aeroelastic computations of flexible configurations, 
Computers &: Structures 30 (1988) 15-28. 

[11] R. D. Rausch, J. T. Batina and T. Y. Yang, Euler flutter analysis of airfoils 
using unstructured dynamic meshes, AIAA Paper No. 89-13834, 30th Struc- 
tures, Structural Dynamics and Materials Conference, Mobile, Alabama, April 
3-5, 1989. 

[12] M. Blair, M. H. Williams and T. A. Weisshaax, Time domain simulations of 
a flexible wing in subsonic compressible flow, AIAA Paper No. 90-1153, AIAA 
8th Applied Aerodynamics Conference, Portland, Oregon, August 20-22, 1990. 

[13] T. W. Strganac and D. T. Mook, Numerical model of unsteady subsonic 
aeroelastic behavior, AIAA Journal 28 (1990) 903-909. 

[14] E. Pramono and S. K. Weeratunga, Aeroelastic computations for wings 
through direct coupling on distributed-memory MIMD parallel computers, AIAA 
Paper No. 94-0095, 32nd Aerospace Sciences Meeting &: Exhibit, Reno, January 
10-13, 1994. 

[15] T. J. R. Hugh es and W. K. Liu, Implicit-explicit finite elements in transient 
analysis: stability theory, ASME Journal of Applied Mechanics 45 (1978) 371- 
374. 

[16] T. Belytschko, P. Smolenski and W. K. Liu, Stability of multi-time step 
partitioned integrators for first-order finite element systems, Comput. Meths. 
Appl. Mech. Engrg. 49 (1985) 281-297. 

[17] Y. C. Fung, An introduction to the theory of aeroelasticity, Dover Publica- 
tions, Inc., New York (1969) 160-185. 

[18] H. Schlichting, Boundary layer theory, Fourth edition, McGraw-Hill, New 
York, 1960. 

[19] T. J. R. Hughes, Analysis of transient algorithms with particular reference 
to stability behavior, T. Belytschko and T. J. R. Hughes, Eds., North-Holland 
Pub. Co. (1983) 67-155. 

[20] C. Farhat, K. C. Park and Y. D. Pelerin, An unconditionally stable staggered 
algorithm for transient finite element analysis of coupled thermoelastic problems, 
Comput. Meths. Appl. Mech. Engrg. 85 (1991) 349-365. 

[21] C. Farhat and N. Sobh, A consistency analysis of a class of concurrent tran- 
sient implicit/explicit algorithms, Comput. Meths. Appl. Mech. Engrg. 84 
(1990) 147-162. 

[22] C. Farhat, S. Lanteri, N. Maman, and P.S. Chen, Parallel and heterogeneous 
solution of three-dimensional transient coupled aeroelastic problems,” (submitted 
for publication) 


51 



[23] N. Maman and C. Faxhat, Matching fluid and structure meshes for aeroelastic 

computations: a parallel approach, Report CU-CSSC-93-12, The University of 
Colorado at Boulder, June 1993. ... 

[24] L. Crivelli and C. Farhat, Implicit transient finite element structural compu- 
tations on MIMD systems: FETI v.s. direct solvers, AIAA Paper No. 93-1310, 
AIAA 34th Structural Dynamics Meeting, La Jolla, California, April 19-21, 1993. 

[25] C. Farhat, L. Crivelli and F. X. Roux, Extending substructure based iterative 
solvers to multiple load and repeated analyses, Comput. Meths. Appl. Mech. 
Engrg. (in press) 

[26] C. Faxhat, L. Fezoui, and S. Lanteri, Two-dimensional viscous flow compu- 
tations on the Connection Machine: unstructured meshes, upwind schemes, and 
massively parallel computations, Comput. Meths. Appl. Mech. Engrg. 102 
(1991) 61-88. 

[27] E. Baxszcz, Intercube communication on the iPSC/860, Scalable High Per- 
formance Computing Conference, Williamsburg, April 26-29, 1992. 

[28] S. Piperno, Numerical methods used in aeroelasticity simulations, Rapport 
CERMICS 92-5, Oct. 1992. 


APPENDIX II 


Simulation of Compressible Viscous Flows on a 
Variety of MPPs: Computational Algorithms for 
Unstructured Dynamic Meshes and Performance Results 


by 


C. Farhat and S. Lanteri 


Department of Aerospace Engineering Sciences and 
Center for Space Structures and Controls 
University of Colorado 
Boulder , Colorado 80309-0429 


January 1994 


Abstract 

Here we report on our effort in simulating unsteady viscous flows on the iPSC-860, the 
CM-5, and the KSR-1 MPPs (Massively Parallel Processors), using a Monotonic Upwind 
Scheme for Conservation Laws finite volume/finite element method on fully unstructured 
fixed and moving grids. We advocate mesh partitioning with message passing as a portable 
paradigm for parallel processing. We present and discuss several performance results ob- 
tained on all three MPP systems in terms of interprocessor communication costs, I/O, 
scalability, and sheer performance. 
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1 Introduction 

In this paper, we detail our approach to the simulation of large scale, 
steady and unsteady, compressible viscous flows on massively parallel 
processors. We consider the numerical solution of the two-dimensional 
Navier-Stokes equations using a mixed finite element/finite volume 
formulation based on unstructured triangular meshes. The spatial 
approximation method combines a Galerkin centered approximation 
for the viscous terms, and a Roe upwind scheme for the computation 
of the convective fluxes. Higher order accuracy is achieved through 
the use of a piecewise linear interpolation method that follows the 
principle of the MUSCL (Monotonic Upwind Scheme for Conservative 
Laws) procedure. The temporal solution is carried out via a 3-step 
variant of the explicit Runge-Kutta method which lends itself to paral- 
lel processing. An ALE (Arbitrary Lagrangian Eulerian) formulation 
is incorporated in the fluid solver to allow the grid points to displace 
in Lagrangian fashion, or be held fixed in Eulerian manner, or be 
moved in some specified way to give a continuous and automatic re- 
zoning capability, depending on the needs of the physical problem to 
be solved. 

Explicit solvers are naturally amenable to parallel processing be- 
cause they essentially involve local computations on vertices, and/or 
edges, and/or triangles of a mesh. However, unstructured meshes in- 
duce indirect addressing memory operations that are costly on many 
hardware architectures. In particular, the present mixed finite ele- 
ment/finite volume solver incurs multiple gather/scatter operations 
between vertex and triangle based arrays. Therefore, in this paper we 
highlight the impact of irregular data access patterns on the compu- 
tational scalability of a parallel unstructured solver, and emphasize 
the importance of data locality in achieving high level performances. 
These concerns and our drive for developing a portable code have led 
us to adopt mesh partitioning with message-passing as a paradigm for 
parallel processing. 

The remainder of this paper is organized as follows. Section 2 
describes the mathematical model of the problem and the approxima- 
tion methods involved in the numerical solution algorithm. Section 3 
identifies the main computational kernels, and motivates the selected 
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parallelization strategy. Overlapping and non-overlapping mesh parti- 
tions are presented, discussed, and contrasted. Finally, Sections 4 and 
5 report and analyze the performance results obtained on the iPSC- 
860, the KSR-1, and the CM-5 parallel processors for various external 
and internal viscous flow simulations, with fixed and moving meshes. 
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2 Simulation of compressible viscous 
flows 

We are interested in the numerical simulation of two-dimensional com- 
pressible viscous flows around or within, fixed, or moving and deform- 
ing bodies. Here, we overview the spatial and temporal discretization 
methods that have been previously detailed in Farhat, Fezoui and 
Lanteri [3] for fixed meshes, and outline the mesh updating procedure 
adopted for aeroelastic computations. 


2.1 Governing equations 

Let Q. C 1R 2 be the flow domain of interest and T be its boundary. 
The conservative law form of the equations describing two-dimensional 
Navier-Stokes flows is given by : 

+ v.f(w(x,t)) = (i) 

where x and t denote the spatial and temporal variables, and 
W = (p,pu,pv,Ef , V=(|,|) 

and 

F(W) and G ( W ) denote the convective fluxes and are given by: 
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puv 
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while R(W ) and 5(W) denote the diffusive fluxes and are given by: 
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xy 
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2 SIMULATION OF COMPRESSIBLE VISCOUS FLOWS 


In the above expressions, p is the density, U = (u,v) is the velocity 
vector, E is the total energy per unit of volume, p is the pressure, e 
is the specific internal energy, r XXl r ry , and r yy are the components of 
the two-dimensional Cauchy stress tensor, k is the normalized thermal 

conductivity, Re — ^° - where po, Uo, Lo and po denote the 

characteristic density, velocity, length, and diffusivity is the Reynolds 

number, and Pr = is the Prandtl number. 

ko 

The velocity, energy, and pressure are related by the equation of 
state for a perfect gas: . r.- 

p = ( T -l)(£-y II fff) 

where 7 is the ratio of specific heats (7 = 1 A for air), and the specific 
internal energy is related to the temperature via: 

£ = C V T = f - \ || U || 2 
P 2 

The components of the Cauchy stress tensor are related to the veloc- 
ities via: 


2 f n du dv\ 2 ( dv du\ 

- V Vai v % ' fcj ’ 



where fi denotes the normalized viscosity. 


In order to account for a potential motion or deformation of the 
computational grid, the following coordinate transformation is intro- 
duced: 


I 


— * 


= £(*,<) 

= t 


( 2 ) 


Throughout this paper, it is assumed that the Jacobian of the above 
transformation defined as: 


5 
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does not vanish at any point and any time. Introducing the grid 
velocity: 


- dx. 

W = r 

dt '< 

and using Eqs. (2), Eq. (1) can be transformed into the following 
ALE (Arbitrary Lagrangian Eulerian) formulation (see, for example, 
Donea [2]): 


d(JW), 

dt 1 


c + Jv.f c (W) = j^v.n(W) 


( 3 ) 


where : 


and F C (W ) and G C {W ) are the modified convective fluxes given by: 



( pu \ 


( pv ^ 

F C (W) = 

puu + p 

, G C (W) = 

puv 


puv 


pvv + p 


\ Eu + up J 


^ Ev + pv J 

and: 





J u 

= u - w x 



t* 

= V~Wy 



( 4 ) 


2.2 Boundary conditions 

The boundary T(t) of the flow domain is partitioned into a wall bound- 
ary r„,(t) and an infinity boundary ^(t): T(t) = r„,(t) U T^t). Let 
n(t ) denote the outward unit normal at any point of T(t), and U w and 
T w denote the wall velocity and temperature. 

On the wall boundary r„,(t), a no-slip condition and a Dirichlet 
condition on the temperature are imposed: 

U = U W , T = T W (5) 

No boundary condition is specified fpr the density. Hence, the total 
energy per unit of volume and the pressure on the wall are given by : 
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p=(l-l)pC v T w , E = pC v T w + \ P \\U W || 2 (6) 

For external flows around airfoils, the viscous effects are assumed 
to be negligible at infinity, so that a uniform free-stream state vector 
Woo is imposed on T^t): 



where a is the angle of attack, and Moo is the free-stream Mach num- 
ber. 

For internal flows, r co (f) is partitioned into upstream and down- 
stream boundaries which are in general in contact with the wall bound- 
ary: T c>0 (<) = rs(t) U In that case, the previous definition 

of Woo is improved using a parabolic profile for the free-stream veloc- 
ity. For example, for the horizontal flow between two plates shown in 
Figure 1, one can specify: 

Uoo = (o , v 00 (y)) T (8) 



Figure 1: Horizontal velocity profile for internal flows 


2.3 Spatial discretization 

The flow domain H(t) is assumed to be a polygonal bounded region 
of ffi, 2 . Let T h be a standard triangulation of Q.(t), and h the maximal 
length of the edges of Th. A vertex of a triangle T is denoted by 5,-, 
and the set of its neighboring vertices by K{I). At each vertex 5,, 
a cell Ci(t) is constructed as the union of the subtri angles resulting 
from the subdivision by means of the medians of each triangle of Th 
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that is connected to Si (see Figure 2). The boundary of Ci(t) is 
denoted by dCi(t ), and the unit vector of the outward normal to dCi(t) 
by #(t) = (j'jx(f)j ^‘y(0)* The union of all these control volumes 
constitutes a discretization of domain £l(t) : 

ns 

«»(<) = U 

t=l 



Figure 2: Control volume in an unstructured grid 

The spatial discretization method adopted here combines the fol- 
lowing features (see [3] for details): 

• a finite volume upwind approximation method for the convec- 
tive fluxes. Second order spatial accuracy is achieved using an 
extension of Van Leer’s MUSCL technique [13] to unstructured 
meshes; 

t a classical Galerkin finite element centered approximation for the 
diffusive fluxes. 

Let Cf and Cf denote the two representations of C,(f) in the co- 
ordinate systems defined by x and <f, respectively. By definition, Cf 
is a fixed reference representation of C,-(t). Integrating Eq. (3) over 

Cf yields: 
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/ / + / / JiJiyrvz- S / m 

cf cf cf 

Given that the time derivative is computed for a constant f and that 
Cf does not depend on the time t, the mapping = £*(x, J) and the 
identity dx = can be used to transform the left-hand side of Eq. 
(9) into: 


1 / }**+! f±*.nyr) m 

Cf cf cf 

Finally, integrating Eq. (10) by parts leads to: 


Ttll Wd • + 

Ci(t) 

£ J T -(W).Vida 

i€K(i) 9C f } (t) 

< 1 > 

+ 

J F c {W).nido 

< 2 > 


ac,(t)nr w (t) 


+ 

J f c {W).n t da 

< 3 > 

= 

9C,(<) nr»(t) 

j, 

< 4 > 


where dCfj(t) = dCf(t)DdCf(t ), and N? = N?(x, y) is the PI shape 
function defined at the vertex 5,* and associated with the triangle T. 

A first order finite volume discretization of < 1 > goes as follows: 


< 1 > = v4" +1 W l n+1 “ A?W? + A t Y, *X, (12) 

j€ K(i) 

where .4’“ denotes the area of the control volume Cf measured at 
time t n , Pij(t) denotes a spatial mean value the normal to dCij(t) y the 
tilde superscript designates a temporal mean value between t n and 
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Z n+1 , and denotes a numerical flux function that approximates 
the following quantity: 

t n + 1 

A ttMWuWjMi)** j J ? e ( W)M { da (13) 

t n dCf } 

Upwinding is introduced by extending Roe’s approximate Riemann 
solver [11] to dynamic meshes and computing $jr c as follows: 


*F'{W»W h Pij) 


2 M ' 3 

- | A(W it Wj, Vij) - (w.Vij)I | — ~ — 


(14) 


where A is Roe’s mean value of the flux Jacobian matrix 


df(W) 

dW ’ 


I is the identity matrix, and the dot product w.Vij is computed as 
suggested recently by N’Konga and Guillard [10]: 


w.Vij = 2 + v(P2ij)) Aj 


(15) 


where Puj and Puj are the end points of the bi-segment dC tJ (see 3). 



Figure 3: Computation of w.Uij 

The control volume area .4” +1 is updated using the following finite 
volume scheme: 


A[* +1 = A? + A*, £ tiJ tj 
}SK(i) 


( 16 ) 



10 
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Following the MUSCL technique, second order accuracy is achieved 
in expression (13) via a piecewise linear interpolation of the states Wij 
and Wji at the interface between cells C, and C } . This requires the 
evaluation of the gradient of the solution at each vertex as follows: 


' W* =w? + i(Viy)f .5^- 



W* = w; - ^w^.srs, 


(17) 


where W* = (p , u , v , p) T — in other words, the interpolation is 
performed on the physical variables instead of the conservative vari- 
ables. The approximate nodal gradients (^WO^ are obtained using 
a /3-combination of centered and fully upwind gradients : 


(VkF)f = (1 - (3)(VW)? ent + PiVW)? 1 ™ (18) 

The half-upwind scheme (/3 = i) is simply obtained by means of a lin- 
ear interpolation of the Galerkin gradients computed on each triangle 
ofC,: 



and the centered 
that verifies: 


J f VW\ T dx 

c, 

1 yareaiT) X W *$ N T 

area(Ci ) ^ 7 - 


(19) 


TiC{ “ Jk=l,Jfc«T 
gradient (VW')f ent (/3 = 0) is given by any vector 


(VVF) f ent .S~Sj - Wj - Wi 


( 20 ) 


The second term < 2 > and the third term < 3 > of Eq. (11) 
include the contributions of the boundary conditions and are evaluated 
as follows: 
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Wall boundary : the no-slip condition is enforced with a strong 
formulation and therefore the corresponding boundary integral in < 
2 > is not explicitly computed. 

Inflow and outflow boundaries : at these boundaries, a precise set 
of compatible exterior data that depend on the flow regime and the 
velocity direction must be specified. Here, a plus-minus flux splitting 
is applied between exterior data and interior values. More specifically, 
the boundary integral < 3 > is evaluated using a non-reflective version 
of the flux-splitting of Steger and Warming [12]: 

J F{W).n t da = A + (Wi,i7 ioo ).W t + (21) 

ac,(t)nroo 


Finally, the viscous integral < 4 > is evaluated using a classical 
Galerkin finite element PI method. The components of the stress 
tensor and those of VN? are constant in each triangle. The velocity 
vector in a triangle is computed as follows: 

Ut = \ £ V k 

4 k—l keT 

Consequently, the viscous fluxes are approximated as follows: 


£ [ n{W).VNfdx = £ area(T)(R T ^f + S T ^p\ 

T,SieTj- T,SieT \ ox oy j 

where Rj and St are the constant values of R{W) and S(W) in the 
triangle T. 


2.4 Time integration 

The resulting semi- discrete fluid flow equations can be written as: 

dW . . 

— + i>(W) = 0 

Because it lends itself to massive parallelism, the following 3-step 
variant of the explicit Runge-Kutta algorithm is selected for time in- 
tegrating the above equations: 
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' w(o) 

< wW 

( 22 ) 

This scheme is often referred to as a low-storage Runge-Kutta algo- 
rithm because only the solution at substep A: — 1 is needed to compute 
the one at substep k. The coefficients are the standard Runge- 
Kutta coefficients and are given by: 

1 


The above time integration algorithm is third order accurate in the 
linear case, and second order accurate in the general non-linear case. 


= W n = W(t = nAt) 

= *= 1 . 2,3 

— w n + 1 


2.5 Dynamic meshes 

In this work, an unstructured dynamic fluid mesh is represented by a 
pseudo structural model (see, for example, Batina [1]) where a ficti- 
tious linear spring is associated with each edge connecting two fluid 
grid points Si and Sj and is attributed the following stiffness: 

kij = i = ( 23 ) 

ySi " x *) 2 + (yj - Vi ) 2 

The grid points located on the downstream and upstream boundaries 
are held fixed. The motion of those points located on the wall bound- 
ary is determined from the wall motion and/or deformation. At each 
time step the new position of the interior grid points is deter- 

mined from the solution of a displacement driven pseudo structural 
problem via a two-step iterative procedure. First, the displacements 
of the interior grid points are predicted by extrapolating the previous 
displacements at time steps t n and t n ~ l in the following manner: 

f hi = 6 n - l Zi 


hi = 2 6 n yi - 8 n - l Vi 


( 24 ) 



mi il Lilli 
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with 8 n x = x n+1 — x n . Next, the above predictions are corrected with 
a few explicit Jacobi relaxations on the static equilibrium equations 
as follows: 


S n+1 Xi 


E kijSij 

E k 'i 

i€*(0 _ 

E kjhj 


s n+l i n 


i6*(0 

E *«; 


l JeK(i) 

Finally, the new positions are computed as: 


x” +1 =X? + ^ n+1 I,- 

y? +1 = y? + ^ n+1 y« 


(25) 


(26) 


U 
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3 Computational and parallel imple- 
mentation issues 

3.1 Identification of the computational ker- 
nels 

From Eqs. (11-22), it follows that our fluid solver contains essentially 
two kernels of elementary computations, one for the convective fluxes, 
and the other for the diffusive ones. Both type of computations can be 
described as three-step sequences of the form Gather-Compute-Scat'ter, 


3.1.1 The convective flux kernel 

The evaluation of the second term of < 1 > in Eq. (11) using the 
numerical flux function $;r c (14) with the second order approximation 
outlined in Eq. (17) can be summarized as follows: 


' II X1 = J ? c (W).v i da = $r c , J (W„W j ,v ij ) 

* dC,j(t) 

k Hji = Rij 

where: 


(27) 


Vi ,j(t) = J Vida - i7i (t) + v 2 (t) 


(28) 


Essentially, one-dimensional elementary convective fluxes are com- 
puted at the intersection between the control volumes C,(f) and Cj(t) 
(see Figure 4 below). Each elementary flux contributes to a flux bal- 
ance at the boundary of the control volume Ci(t). This balance in- 
volves the accumulation over the set of neighboring vertices K(i) of 
all computed fluxes. From the second of Eqs. (27), it follows that 
only Hij needs to be computed in order to update the flux balances 
at the two end-point nodal values of edge Eij = {Si, Sj} . Therefore, 
the most efficient way for evaluating the convective fluxes is to loop 
over the list of the mesh edges and compute as follows: 
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Figure 4: Evaluation of a convective flux along an edge {5 1 ,, Fj} 


R1 

Irk 


For each edge Eij = of 7^ Do 

Gather W, = W(S,-) , W { = TF(S,) 

Gather VW t = VW(Si) , VW; = W(5.) 

Compute Hij 
Scatter $,• = + H;j 

Scatter — H ,j 

End Do 

3.1.2 The diffusive flux and nodal gradient kernel 

In the last term < 4 > of Eq. (11), the elementary diffusive flux 
TZi(T) is constant in each triangle T. Its evaluation requires accessing 
the values of the physical state W at the three vertices 5, , Sj and 


f = J 1 1i(W).VNT dxdy 


ON? on t 

= area(T)(R T ^- + S T ^-) 


(29) 


dx dy 

The values of R? and St contribute tp the diffusive fluxes at all three 
vertices of triangle T. The sum symbol in < 4 > is a clear indication 
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of a gather operation. Clearly, the most efficient way for evaluating 
the convective fluxes is to loop over the list of the mesh triangles and 
compute as follows: 

For each element T{j k = {•5'i, Sj, 5^} of Do 

Gather W t = W(Si) , Wj = VF(Sj) , W k = 

W (S k ) 

Compute Rt, St 
Scatter V,- = V,- + Hi{T) 

Scatter Vj = Vj + Hj{T) 

Scatter V k = V k + Rk(T) 

End Do 

The evaluation of the half- upwind nodal gradient (19) follows the 
same computational pattern described above. 

3.2 The mesh partitioning with message-passing 
parallel paradigm 

In addition to efficiency and parallel scalability, portability should be 
a major concern. With the proliferation of computer architectures, 
it is essential to adopt a programming model that does not require 
rewriting thousands of lines of code — or even worse, altering the 
architectural foundations of a code — every time a new parallel pro- 
cessor emerges. Here, we are neither referring to differences between 
programming languages, nor to differences between the multitude of 
parallel extensions to a specific programming language. We are more 
concerned about the impact of a given parallel hardware architecture 
on the software design, and sometimes, on the solution algorithm it- 
self. For example, a data parallel code written for the CM-2 or CM-5 
machines could require major rehauling before it can be adapted to an 
iPSC computer. A parallel-do-loop based code can be easily ported 
across different true shared memory multiprocessors, but may require 
substantial modifications before it can run successfully on some dis- 
tributed memory systems. 

Based on our w hands on” experience with a dozen of different par- 
allel processors, we believe that the fnesh partitioning and message- 
passing lead to portable software designs for parallel computational 


3 COMPUTATIONAL AND PARALLEL IMPLEMENTATION ISSUES 1 7 


mechanics. Essentially, the underlying mesh is assumed to be parti- 
tioned into several submeshes, each defining a subdomain. The same 
“old” serial code can be executed within every subdomain. The as- 
sembly of the subdomain results can be implemented in a separate 
software module and optimized for a given machine. This approach 
enforces data locality, and therefore is suitable for all parallel hardware 
architectures. For example, we have shown in [8] that for unstructured 
meshes, this approach produces substantially better performance re- 
sults on the KSR-1 than the acclaimed virtual shared memory pro- 
gramming model. Note that in this context, message-passing refers 
to the assembly phase of the subdomain results. However, it does 
not imply that messages have to be explicitly exchanged between the 
subdomains. For example, message-passing can be implemented on a 
shared memory multiprocessor as a simple access to a shared buffer, 
or as a duplication of one buffer into another one. 

In this work, we use essentially the same code on the iPSC-860, 
the KSR-1, and the CM-5 parallel processors. This code also runs 
on a workstation. We consider mesh partitions with and without 
overlapping for reasons that are discussed next. 

3.2.1 Overlapping mesh partitions 

The reader can verify that for the computations described herein, 
mesh partitions with overlapping simplify the programming of the sub- 
domain interfacing module. Only one communication step is required, 
after the local physical states have been updated. Depending on the 
order of the spatial approximation, the overlapping region can be one 
or three triangles wide (see Figure 5 below). For a first order spatial 
approximation, we have by definition: 


w* 

rr *J 

II 

J* 

r 

II 


which shows that the overlapping region needs in that case to be only 
one triangle wide. 

However, mesh partitions with overlapping also have a drawback: 
they incur redundant floating-point operations. 
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Figure 5: Overlapping mesh partition 
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For fixed meshes and overlapping partitions, the main loop of the 
parallel fluid solver described herein goes as follows: 

Repeat step = step + 1 

Compute the local time steps 
For srk = 1 to nsrk Do 

Compute the nodal gradients 
Compute the diffusive fluxes 
Compute the convective fluxes 
Update the physicalstates 
Exchange the conservatives variables 

End Do 

Until step = step max 

In the above pseudo code, step max denotes the maximum number of 
time steps, and nsrk denotes the number of steps in the Runge-Kutta 
integration algorithm. 

All overlapping mesh partitions used in this investigation were gen- 
erated by the decomposer described in [9]. 

3,2.2 Non-overlapping mesh partitions 

Non-overlapping mesh partitions (see Figure 6 below) incur little re- 
dundant floating-point operations but induce one additional communi- 
cation step. While physical state variables are exchanged between the 
subdomains in overlapping mesh partitions, partially gathered nodal 
gradients and partially gathered fluxes are exchanged between subdo- 
mains in non-overlapping ones. 
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Figure 6: Non-overlapping mesh partition 

For fixed meshes and non-overlapping partitions, the main loop of 
the parallel fluid solver described herein goes as follows: 

Repeat step = step + 1 

Compute the local time steps 
For srk = 1 to nsrk Do 

Compute the nodal gradients 
Compute the diffusive fluxes 
Exchange the nodal gradients 
Compute the convective fluxes 
Exchange the convective fluxes 
Update the physical states 

End Do 

Until step = step max 
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All non-overlapping mesh partitions discussed in this investigation 
were generated by the TOP/DOMDEC software described in [4], 

3.2.3 To overlap or not to overlap? 

To answer this question, we analyze the communication requirements 
of both families of mesh partitions, and the amount of redundant com- 
putations they incur for two-dimensional fixed problems. Because we 
are interested in a comparative study, it suffices to consider the case 
of a single interface between two subdomains with uniform triangula- 
tions. 

Let denote the number of interface vertices in a non overlap- 
ping mesh partition (see Figure 7). In the first communication step, 
8 X n 1 j ov words related to the nodal gradients are exchanged between 
the two subdomains. In the second communication step, 4xn” ov fluxes 
are exchanged. Hence, the total communication cost per subdomain 
is given by: 


T™ = 2xT s + l2x nY v x T t (31) 

where T s denotes the startup time of a message, and T r denotes the 
transmit time for a 64-bit word. 

In a non-overlapping mesh partition, the only redundant computa- 
tions are those associated with the evaluation of the convective fluxes 
along the interface edges. Since an elementary convective flux requires 
about 200 floating-point operations, the total time per subdomain as- 
sociated with redundant computations can be estimated as: 

T™£ = 200 x {nf v - 1) x T« (32) 

where T a denotes the times it takes to perform a single floating-point 
operation. 

From Figure 8, it follows that the total number of overlapping 
vertices in an overlapping mesh partition is 4 x n* ot \ In this case, 
only one communication step is required to exchange 4 components of 
the physical state at half of the overlapping vertices. Hence, the total 
amount of communication per subdoihain is given by: 
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D2 
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O 


vertices updated by D2 and communicated to Dl 


Figure 8: Analysis of an overlapping mesh partition 
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T™ m = T 3 + 8 x rtf'* x T r (33) 

In an overlapping mesh partition, redundant computations are 
performed during the evaluation of both the convective and diffu- 
sive fluxes. Given that an elementary diffusive flux requires about 
100 floating-point operations, the total time per subdomain associ- 
ated with redudant computations is given in that case by: 

Trcd = (200 x (4 x n 1 }™) + 100 x (4 X ( n] ov - 1)) X T a (34) 

where 4 X 7 i* ov is the number of overlapping edges that generate redun- 
dant convective flux computations, and 4 x (nj ou — 1) is the number 
of overlapping triangles that generate redundant diffusive flux compu- 
tations. 

From Eqs. (31,33), it follows that for sufficiently large messages 
we have: 

IZn = r, + 8xnr xr r W 2 

T c n ™ 2 x T a + 12 x n n j ov xT r ~ 3 ^ J 

which shows that overlapping the mesh partitions reduces the com- 
munication costs by 33%. 

However from Eqs. (31-34), it follows that for a sufficiently large 

n nov we h ave; 


(T“„+T;a-(T"“+r”7) » 1000 x nj~ X r „(l - x ) (36) 

which suggests that overlapping the mesh partitions will incur a greater 
total parallel overhead than not overlapping them, because of the 
resulting redundant computations. For example on the iPSC-860, 
T a = 0.2 x 10 -6 seconds (sustained 5 Mflops per processor), T r = 
3.184 x 10 -6 seconds (sustained 2.5 Mbytes/second), and (T c ° u m + 
T™ d ) — (T£,™ + T r "™) » 0.18 x 10 -3 x nj v > 0. Nevertheless, this 
specific result also suggest that for two-dimensional problems, simi- 
lar performance results will be obtained for mesh partitions with or 
without overlapping. 
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Finally, we caution the reader that different conclusions may be 
drawn for three-dimensional problems where overlapping can also re- 
quire significantly more storage. 
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4 Performance results on a variety of 
MPPs 

In this section, we discuss the parallel performance results obtained 
on various configurations of the iPSC-860, the KSR-1, and the CM-5 
parallel processors. 


4.1 Focus problem 

We consider the numerical simulation of the unsteady viscous flow 
around a fixed NACA0012 airfoil, starting impulsively from a uniform 
flow. The angle of attack is set to 30°, and the free stream Mach 
number to 0.1. Several physical solutions of this problem were previ- 
ously reported in [3] for different Reynolds numbers. All performance 
results reported herein are for 100 iterations and 64-bit arithmetic. 
More importantly, the redundant floating-point operations are 
not counted when evaluating the mflop rate, which is a strict ap- 
proach to benchmarking. 

A partial view of an unstructured triangulation of the computa- 
tional domain is given in Figure 9. Seven meshes with increasing sizes 
have been generated. Their characteristics are summarized in Table 
1 below where Nv denotes the number of vertices, Nt the number of 
triangles, and Ne the number of edges. 


MESH 

N v 

N t 

N e 

Wl 

8119 

15998 

24117 

t; : ]M2 

16116 

30936 

47052 

M3 

32236 

63992 

96228 

M4 

63974 

127276 

191250 

M5 

131035 

261126 

392161 

M6 

262717 

523914 

786631 

Ml 

523196 

1044504 

1567700 


Table 1 : seven meshes and their characteristics 


Throughout the remainder of this' paper, the following nomencla- 
ture is used for the investigated mesh partitioning algorithms: 



f 
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Figure 9: Partial view of a NACA0012 mesh 

• SCT : Sector [9] 

• RIB : Recursive inertial bisection [5] 

• GRD : Greedy [5] 

• RGB : Recursive graph bisection [16] 

• RSB : Recursive spectral bisection [16] 

4.2 Parallel scalability for increasing size prob- 
lems 

Parallel scalability is evaluated here for problems where the subdomain 
size is fixed, and the total size is increased with the number of pro- 
cessors. Note that because we are dealing with unstructured meshes, 
some slight deviations are inevitable. Overlapping mesh partitions are 
generated using the RIB heuristic. 

Tables 2-3 summarize the performance results obtained on the 
iPSC-860 and KSR-1 parallel systems. The number of processors is 
denoted as N p . The parallel CPU time and the Mflop rate are shown 
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to remain almost constant when the problem size is increased with 
the number of processors, which demonstrates the scalability of the 
parallel solver. The slight degradations in efficiency are mainly at- 
tributed to overlapping since redundant operations are not accounted 
for in the evaluation of the mflop rate. The KSR-1 processor cell is a 
RISC-style superscalar 64-bit unit operating at a peak of 40 Mflops. 
Clearly, despite a rather large number of gather/scatter operations, 
25% of this peak performance Is attained, — - 


N v 

N ? 

CPU Time 

Mflop/s 

Comm Time 

% Comm 

8119 

1 

491.2 s 

6 

0 s 

0 

16116 

2 

491.7 s 

11 

28.6 s 

5.81 

32236 

4 

514.7 s 

22 

31.3 s 

6.09 

63974 

8 

519.5 s 

43 

40.0 s 

7.71 

131035 

16 

536.6 s 

85 

40.3 s 

7.52 

262717 

32 

548.6 s 

159 

44.8 s 

8.17 


Table 2 : Parallel scalability for increasing size problems 
Computations with overlapping mesh partitions on the iPSC-860 


N v 

N p 

CPU Time 

Mflop/s 

Comm Time 

% Comm 

8119 

1 

272.3 s 

10 

0 s 

0 

16116 

2 

272.3 s 

20 

1.8 s 

0.67 

32236 

4 

286.0 s 

39 

3.2 s 

1.18 

63974 

8 

289.3 s 

77 

4.3 s 

1.52 

131035 

16 

306.7 s 

149 

5.9 s 

1.95 

262717 

32 

316.1 s 

276 

8.4 s 

2.66 


Table 3 : Parallel scalability for increasing size problems 
Computations with overlapped mesh partitions on the KSR-1 


4.3 Influence of the mesh partitioning algo- 
rithm 

Next, we focus on mesh M6 with 32 processors and non-overlapping 
partitions, and investigate the influence of the partitioning algorithm 
on parallel performance. Tables 4 and 5 report the measured CPU 
time and Mflop rates. “Conv” and “DifF designate respectively the 
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convective and diffusive fluxes. In all cases, the RSB algorithm yields 
the fastest solution time, even when it does not produce the smallest 
communication time. The reason is that, for mesh M6, the RSB al- 
gorithm does a better job than the others at generating subdomains 
that are well balanced vertex-wise, element-wise, and edge-wise, si- 
multaneously. 


Decomp 

CPU Time 

Conv Time 

Diff Time 

Comm Time 

Mflop/s 

GRD 

550.0 s 

330.4 s 

163.4 s 

47.4 s 

158 

RGB 

556.3 s 

335.1 s 

157.9 s 

62.7 s 

158 

RSB 

538.0 s 

326.9 s 

153.7 s 

53.8 s 

162 


Table 4 ^Influence of the mesh partitioning algorithm 
Computations with non-overlapping mesh partitions on an iPSC-860/32 


Decomp 

CPU Time 

Conv Time 

Diff Time 

Comm Time 

Mflop/s 

GRD 

343.7 s 

165.9 s 

112.8 s 

18.9 s 

254 

RGB 

340.1 s 

168.3 s 

105.1 s 

18.5 s 

256 

RSB 

322.5 s 

160.4 s 

101.4 s 

14.2 s 

270 


Table 5 : Influence of the mesh partitioning algorithm 
Computations with non-overlapping mesh partitions on a KSR-1/32 


The reader can verify that the numbers reported in columns 3 to 5 
of Tables 4-5 do not add up to the total CPU time reported in column 
2. The difference corresponds to various parallel and sequential over- 
heads. On the iPSC-860, these overheads represent less than 2% of the 
total solution time. However on the KSR-1, they represent about 14% 
of the total CPU time. Indeed, the computing mode on the iPSC-860 
is parallel by default, while on the KSR-1 it is sequential by default. 
Hence, many fork- join type of procedures are necessary on the KSR-1, 
which explains the relatively large amount of overhead. 

Also, note that using the same number of processors, the KSR-1 is 
twice as fast as the iPSC-860 at computing the convective fluxes, but 
only 1.5 times faster at computing the diffusive ones. This is because 
the evaluation of the convective fluxes requires less indirect addressing 
than the evaluation of the diffusive ohes. 
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Figure 12: Vertex-wise load bal. - 32 non-overlapping sub. 

RSB algorithm 
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Figure 13: Element-wise load bal. - 32 non- overlapping sub. 

GRD algorithm 
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Figure 16: Edge-wise load bal. - 32 non-overlapping sub. 
GRD algorithm 



Figure 17: Edge- wise load bal. - 32 non-overlapping sub. 
RGB algorithm 





4 PERFORMANCE RESULTS ON A VARIETY OF MPPS 


34 



Figure 18: Edge-wise load bal. - 32 non-overlapping sub. 

RSB algorithm 


4.4 Parallel scalability for fixed size problems 

Tables 6-8 summarize the performance results obtained on the iPSC- 
860 and KSR-1 parallel processors for overlapping mesh partitions 
generated by the RIB algorithm. For the largest mesh M 7, a Gi- 
gaflop performance level is attained using 128 processors of the KSR-1 
system. Good scalability is observed on both machines. 


N ? 

CPU Time 

Conv Time 

Diff Time 

Comm Time 

Mflop/s 

32 

548.7 s 

323.1 s 

164.6 s 

44.9 s 

159 

64 

282.7 s 

163.6 s 

82.4 s 

24.4 s 

309 

128 

144.5 s 

82.2 s 

42.7 s 

18.9 s 

617 


Table 6 : Performance results for mesh M6 
Computations with overlapping mesh partitions on the iPSC-860 
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N f 

CPU Time 

Conv Time 

DifF Time 

Comm Time 

Mflop/s 

32 

316.1 s 

155.4 s 

111.3 s 

8.4 s 

276 

64 

169.7 s 

79.0 s 

58.1 s 

9.8 s 

514 

128 

93.0 s 

39.7 s 

31.2 s 

6.4 s 

938 


Table 7 : Performance results for mesh M6 
Computations with overlapping mesh partitions on the KSR-1 


N P 

CPU Time 

Conv Time 

Diff Time 

Comm Time 

Mflop/s 

128 

174.0 s 

77.7 s 

60.2 s 

14.0 s 

1024 


Table 8 : Performance results for mesh M 7 
Computations with overlapping mesh partitions on the KSR-1 


Finally, Tables 9-10 compare the performances of the iPSC-860 and 
KSR-1 parallel systems for the case of mesh M6 and non-overlapping 
mesh partitions generated by the GRD algorithm. For the same 
number of processors, the KSR-1 machine is reported to be 1.67 times 
faster than the iPSC-860, even though its basic processor is supposed 
to be 1.5 times slower than that of the iPSC-860. 


| N p 

CPU Time 

Conv Time 

Diff Time 

Comm Time 

Mflop/s 

1 64 

287.7 s 

164.0 s 

79.8 s 

42.3 s 303 


Table 9 : Performance results for mesh M6 
Computations with non-overlapping mesh partitions on the iPSC-860 


N p 

CPU Time 

Conv Time 

Diff Time 

Comm Time 

Mflop/s 

64 

171.8 s 

83.1 s 

51.1 s 

14.1 s 

508 

128 

94.8 s 

40.1 s 

24.3 s 

15.2 s 

921 


Table 10 : Performance results for mesh M6 
Computations with non- overlapping mesh partitions on the KSR-1 


4.5 Performance results on the CM-5 

Recently, we have implemented our fluid solver on a 32 processor CM- 
5 system using Fortran 77 on a node and the CMMD message pass- 
ing library (this corresponds to the Sparc model of computation on 
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the CM-5). For mesh M6 and a 32 subdomain decomposition with 
overlapping using the RIB algorithm, Table 11 reports the measured 
performance results and compares them to those obtained on an iPSC- 
860/32 and a KSR-1/32 computers. 


MPP 


Conv Time 

Diff Time 

Comm Time 

Mflop/s 

CM-5 

855.4 s 

541.3 s 

241.7 s 

44.0 s 

102 

iPSC-860 

548.6 s 

323.1 s 

164.6 s 

44.8 s 

159 

KSR-1 

316.1 s 

155.4 s 

111.3 s 

8.4 s 

276 


Table 1 1 : Performance results for mesh M6 
Computations with overlapping mesh partitions 


Clearly, the results reported in Table 11 for the CM-5 are not as 
impressive as those reported, for example, in [15, ?]. This can be 
attributed to several factors including the use of the message-passing 
model for portability reasons, the variety of gather/scatter operations 
required by our specific fluid solver, and more importantly our strict 
approach to performance benchmarking. 

Next, we report performance results on the CM-5 using a global 
CM Fortran approach (this corresponds to the Vector Units model 
of computation on the CM-5), and the parallel version of our solver 
that was previously developed for the CM-2/200 and described in 
[3]. In this so called data parallel approach, all local computations 
are carried out on a control volume. Therefore, all data structures 
are vertex based. Note that this approach generates a substantial 
amount of redundant computations. The reportecf comm un ication 
timings correspond to the inter and intra vector units gather/scatter 
operations. In Table 12, MFU is a Mflop rate that does not account 
for redundant arithmetic operations, while MFR is a Mflop rate that 
does. 


CPU Time 

Conv Time 

Diff Time 

Comm Time 

MFU 

MFR 

342.0 s 

116.9 s 

49.0 s 

162.0 s 

107 

238 


Table 12 : Performance results for mesh M5 
Computations with overlapping mesh partitions on the CM-5 


























4 PERFORMANCE RESULTS ON A VARIETY OF MPPS 


The results reported in Table 12 show that about 50% of the 
elapsed time is spent in gather /scatter operations, which is consistent 
with the results obtained by other investigators for two-dimensional 
finite element fluid problems [7]. Note that this percentage is higher 
than those observed on the iPSC-860 and KSR-1 systems. 
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5 Applications 

5.1 Steady viscous flow inside a model jet en- 
gine • -• -r ' •- -- - - i : » ■ v 

First, we consider the numerical simulation of a steady viscous flow 
inside a model jet engine. The free stream Mach number and the 
Reynolds number are set respectively to 0.2 and 2000. The com- 
putational grid is illustrated in Figure 19. Its characteristics are 
N v = 12233 , N t = 22936, and N E = 35170. 

This simulation is carried out on the KSR-1 using overlapping mesh 
partitions generated by the RIB algorithm. Here, the pseudo time 
integration is carried out at CFL=1.9 via a four step Runge-Kutta 
method with a?i = 0.11 , <22 = 0.2766 , <23 = 0.5 and a 4 = 1.0. A local 
time step strategy is introduced in order to accelerate convergence. 
After 1527 iterations, the initial residual is reduced by a factor of 10 4 . 
The resulting steady mach lines are depicted in Figure 20. 

Table 13 reports the performance results obtained on the KSR-1, 
and Table 14 summarizes the characteristics of the generated mesh 
partitions. S(p ) denotes the speed-up using p processors. Nj de- 
notes the total number of interface vertices, and Max Nj denotes the 
maximum number of interface vertices per subdomain. The jump of 
communication costs between the case with N p = 2 and that with 
Np = 4 can be attributed to the accidental increase in Max Nj , which 
implies an increase in the maximum message length. The slightly 
superlinear speed-up observed for N p = 4 is not uncommon on the 
KSR-1. If the problem to be solved does not fit exactly into a single 
processor’s memory, the solution on one processor still involves some 
interprocessor communication that is difficult to time by the user. 


N p 

CPU Time 

Mflop/s 

Comm Time 

% Comm 

s(p) 

1 

8946 s 

9 

— 

— 

1.0 

2 

4532 s 

17 

28.7 s 

0.63 

1.9 

4 

2136 s 

37 

49.9 s 

2.33 

4.1 

8 

1168 s 

68 

46.3 s 

3.96 

7.6 


Table 13 : Performance results for an internal flow simulation 
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Computations with overlapping mesh partitions on the KSR-1 


N ? 

Min Nv 

Max Ny 

Min Nt 

Max Nt 

Max Ni 

2 

6222 

6229 

11568 

11668 

113 

4 

3205 

3334 

5942 

6183 

171 

8 

1593 

1774 

2964 

3204 

86 


Table 14 : Characteristics of the mesh partitions 


5.2 Airfoil flutter and control surface 

Next, we consider the flutter simulation of a NACA0012 airfoil in tran- 
sonic flow. The structural dynamics behavior of the airfoil is repre- 
sented by a two-spring two-degree of freedom system. The airfoil twist 
0 is associated with a torsional spring and monitors the angle of attack. 
The lateral deflection h is associated with a lineal spring and monitors 
the airfoil bending. The evolution of this system is governed by a set 
of differential equations that can be written in non-dimensional form 
as follows (for example, see [6]): 

r fit + Elfi + 4 ^ M °° d ± , HMk h = 2M ™ C ‘ 

dP 2 dP V*u$ dt V* 2 u a ftp 

S (37) 

x±£h 2 Hi , dd 7 jMl 2Ml 0 C m 

, 2 dt? 4 dt 2 V * dt V* 2 7 r/i 

In Eqs. (37) above, the bar superscript indicates a non-dimensional 
variable, and C\ and C m denote respectively the lift coefficient and the 
torsional moment. These two quantities are related to the generalized 
aerodynamics forces Qh and Qq by: 

f Qh = j/>ooK£(26)C, 

= J/>ooV2(2i 0 2 c m 


Qe 


( 38 ) 
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:h lines : Min = 0.0 , Max = 0.6 , AM = 0.05 
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where 2b is the airfoil chord, V * = —22. j s the normalized velocity, Vor, 

bu>e 

TCI 

is the free stream velocity, and \i = — is the ratio of the airfoil 

irpoob 2 

mass per transversal unit to the free stream density. 


The effect of an additional control surface such as a flap is sim- 
ulated with the superposition of a controlled motion of a fraction of 
the airfoil trailing edge. This secondary motion is defined by the flap 
angle which obeys the following control law: 

6(t) = G h h(t je^n + Gedityv* (39) 

where G & et G$ are gain coefficients, and <ph and phase angles (see 
Figures 21-23). 
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Figure 23: a = 12° , 8 = 12° , G* = 1.0 
(mesh is coarsened for clarity) 

Two aeroelastic simulations with and without the control surface 
are performed using mesh M4 and: 

26 = 1, m = 1, 6 = ^ = 0 

ug = Uh = 100 rad/s 

70 = 1.865, /i = 60, aft = —2, 10 = 1.8 

V* = 0.6, Q* = 0.006, Foo = 30.0m/s 

The flow initial conditions are identified with the steady solution at 
Moo = 0-8 and zero angle of attack. After the steady state is reached, 
a perturbation in the angle of attack A0 = 0.01 radian is introduced, 
which causes the airfoil to vibrate and the flow to become unsteady. 
The governing aeroelastic equations (37) are used to predict the dy- 
namic response of the system. The unsteady fluid flow equations are 
time integrated with a global time step strategy. At each time step, 
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Figure 21: a = 0° , <5 = 0° , G e = 0.0 
(mesh is coarsened for clarity) 



Figure 22: a = 12° , S'= —6° , Gg = 0.5 
(mesh is coarsened for clarity) 
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the dynamic mesh is updated with 8 explicit Jacobi relaxations as de- 
scribed in Section 2.5* All computations are run on both the iPSC-860 
and KSR-1 parallel processors. 


Figures 24 and 25 report the evolution in time of the angle of at- 
tack 6 and the lift coefficient C \ . Clearly, when the control surface is 
enabled with Gh = G$ = 0.75 and <ph = (p$ = 7r, a stable aeroelas- 
tic response is observed. When it is disabled, a flutter instability is 
reached. 



Figure 24: Angle of attack 9 (in °) versus physical time t (in seconds) 
Gh = Go = 0.75 , <ph = <p& = 7T 
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For this application, the performance results obtained on the iPSC- 
860 and KSR-1 parallel processors are summarized in Tables 15-16. 
Simulation time is reported for 100 steps and includes sequential I/O 
costs for saving the computed solutions on disk. 


ESI 

Simulation Time 

Flux Comp Time 

Mesh Update 


16 

372.0 s 

243.2 s 

32.4 s 

33.5 s 

32 

208.0 s 

124.0 s 

17.7 s 

17.1 s 

64 

101.0 s 

64.0 s 

10.1 s 

12.8 s 


Table 15 : 100 steps of an aeroelastic simulation with mesh M4 
Computations with overlapping mesh partitions on the iPSC-860 


E39 

Simulation Time 

Flux Comp Time 

Mesh Update 

Comm Time 

16 

326.0 s 

135.0 s 

35.8 s 

18.8 s 

32 

176.0 s 

69.2 s 

19.1 s 

18.9 s 

64 

99.7 s 

35.1 s 

11.0 s 

15.8 s 


Table 16 : 100 steps of an aeroelastic simulation with mesh M4 
Computations with overlapping mesh partitions on the KSR-1 


Both parallel processors are shown to deliver good speed-ups. In- 
terprocessor communication Time varies between 9% and 12.7 % on 
the iPSC-860, and 5.7% and 15.8% on the KSR-1. The cost of up- 
dating the dynamic mesh is about 10% of the total cost only. Note 
that for this simulation, the KSR-1 is not reported to be 1.5 times 
faster than the iPSC-860 for the same number of processors, unlike 
in all previous cases. This suggests that I/O on the KSR-1 is more 
expensive than on the iPSC-860. 
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